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Abstract 


An  extended  Kalman  filter  algorithm  which  incorporates 
an  enhanced  correlator/linear  measurement  model  and  a  nonlin¬ 
ear  target  acceleration  dynamics  model  is  described  for  use 
in  a  pointing  and  tracking  system  for  high-energy  ground- 
based  laser  weaponry.  The  measurement  model  used  in  the  fil¬ 
ter  combines  the  computational  benefits  of  a  correlation  algo¬ 
rithm  with  the  statistical  accuracy  available  from  a  Kalman 
filter.  The  dynamics  model,  based  on  a  constant  turn-rate 
target  acceleration  model,  is  deemed  to  be  a  better  represen¬ 
tation  of  the  true  target  dynamics  than  a  linear  first-order 
Gauss-Markov  target  acceleration  model  for  the  cases  of  in¬ 
terest.  Optical  processing  techniques  are  completely  speci¬ 
fied  for  the  correlation  stage  to  perform  the  required  corre¬ 
lation  in  real  time,  and  the  filter  stage  to  perform  the  lin¬ 
ear  algebra  required  for  the  Kalman  filter.  This  extensive 
use  of  optics  allows  the  development  of  two  tracking  algo¬ 
rithms  based  on  the  same  models:  a  FLIR-constrained  tracker 
with  a  30  Hz  frame  rate  and  an  unconstrained  tracker  with  a 
100  Hz  frame  rate  using  real-time  sensors  in  place  of  the 
FLIR. 


I.  Introduction 


1.1  -  Overview 

The  decision  by  the  Department  of  Defense  to  undertake 
massive  development  of  high  energy  laser  weapons  required 
solutions  to  a  number  of  highly  complex  problems  in  a  wide 
range  of  specialty  areas.  While  some  of  these  solutions  are 
fairly  well  developed,  there  remain  major  issues  yet  to  be 
addressed  in  the  pointing  and  tracking  problem  in  order  to 
achieve  a  fully  acceptable  aiming  system. 

To  understand  what  issues  must  be  addressed  requires  a 
review  of  systems  development  to  date.  Since  all  of  the 
unclassified  efforts  which  are  applicable  to  the  system 
development  within  this  thesis  have  been  generated  at  the 
Air  Force  Institute  of  Technology  (AFIT)  under  the  guidance 
of  Dr.  Peter  S.  Maybeck,  this  review  will  be  limited  to 
these  efforts.  The  solutions  obtained  from  this  previous 
work  will  serve  to  highlight  the  advantages  and  limitations 
available  at  the  starting  point  of  this  effort. 

Additionally,  as  most  control  engineers  lack  a  back¬ 
ground  in  the  field  of  optical  processing,  a  brief  review 
of  the  field  with  particular  emphasis  on  algorithms  usable 
in  Kalman  filter  applications  is  included.  While  many  of 
the  recently  developed  optical  techniques  are  highly  complex 
to  people  outside  of  the  field,  they  warrant  investigation 
as  potential  candidates  for  improving  the  overall  filter 
performance.  This  review,  along  with  the  development  of 
the  specific  optical  architecture  in  Chapter  IV,  will 
demonstrate  that  an  optical  solution  includes  a  number  of 
advantages  specifically  applicable  to  the  limitations 
inherent  to  software  versions  of  the  present  Kalman  filter 
solutions . 

The  remainder  of  this  introduction  will  be  built  on  the 
review  of  solutions  which  exhibit  the  greatest  potential 


towards  an  overall  solution  of  the  pointing  and  tracking 
problem.  The  overall  solution  obtained,  using  the  advantages 
of  an  optically-processed  Kalman  filter  to  offset  previously 
obtained  limitations,  will  then  be  presented  as  the  optimal 
solution  within  the  remainder  of  this  thesis. 

Pointing  and  Tracking  Solution  Review 


The  initial  system  solution  (12,16,20)  generated  at 
AFIT  used  a  four-state  simple  extended  Kalman  filter  to  track 
targets  using  outputs  from  a  forward-looking  infrared  (FLIR) 
sensor  as  measurements.  This  technique  allowed  exploitation 
of  knowledge  unused  by  prior  correlation  tracker  develop¬ 
ments:  size,  shape,  and  motion  characteristics  of  the  target 

atmospheric  jitter  spectral  description;  and  background  and 
sensor  noise  characteristics.  Monte  Carlo  performance  analy¬ 
ses  indicated  an  enhanced  ability  of  the  nonadaptive  filter 
to  track  a  realistic  distant  point  source  target  with  an 
error  standard  deviation  of  0.2  picture  elements  (where  pic¬ 
ture  elements  correspond  to  a  field  of  view  of  20  yrad  by 
20  prad) ,  under  expected  tracking  conditions. 

Although  the  simple  extended  Kalman  filter  approach 
achieved  very  accurate  tracking  performance  under  the  nomin¬ 
ally  assumed  conditions,  robustness  studies  portrayed  a  sig¬ 
nificant  degradation  in  performance  when  the  filter's  inter¬ 
nal  model  did  not  depict  the  target's  intensity  profile  or 
motion  characteristics  well  (7,13.14).  Background  noise 
properties  were  shown  to  be  of  secondary  importance  at  expec¬ 
ted  signal- to-noise  ratios.  These  studies  emphasized  the 
need  for  better  models  and  adaptivity  within  the  filter 
structure . 

Further  analysis  resulted  in  incorporation  of  an 
eight-state  filter  using  a  modified  target  dynamics  model  as 
well  as  online  adaptation  to  target  shape  effects,  changing 
target  motion  characteristics,  and  maximum  signal  intensity 
(7,14).  This  was  shown  to  possess  considerable  performance 
potential  for  highly  maneuverable  targets  despite  background 
clutter . 


Improvements  to  the  solutions  generated  by  the  initial 
research  efforts  were  incorporated  into  the  third  major  AFIT 
system  development  (19,23).  This  investigation  was  based  on 
a  tracker  able  to  handle  "multiple  hot  spot"  targets,  in 
which  digital  (or  potentially  optical)  signal  processing  was 
employed  on  the  FLIR  data  to  identify  the  underlying  target 
shape.  This  identified  shape  was  then  used  in  the  measure¬ 
ment  model  portion  of  the  filter  as  it  estimated  target 
position  offset  from  the  center  of  the  field  of  view.  The 
algorithm  then  used  an  extended  Kalman  filter  to  process  the 
raw  intensity  measurements  from  the  FLIR  to  produce  target 
estimates . 

An  alternate  algorithm  also  developed  within  this  re¬ 
search  effort  used  a  linear  Kalman  filter  to  process  the 
position  indications  of  an  enhanced  correlator  in  order  to 
generate  tracking  estimates.  This  enhancement  was  achieved 
by:  (1)  thresholding  to  reject  lower  values  of  the  corre¬ 

lation  function  so  that  a  simple  centroid  calculation 
yielded  a  result  close  to  the  peak  of  the  correlation  func¬ 
tion;  (2)  incorporation  of  the  dynamics  information  from 
the  Kalman  filter;  and  (3),  use  of  the  online-identified 
target  shape  as  a  template  instead  of  merely  using  previous 
frames  of  data.  Both  algorithms,  under  performance  analyses 
using  various  tracking  environment  conditions  and  different 
choices  of  design  parameters,  exhibited  significant  poten¬ 
tial,  with  the  correlation/filter  tracker  involving  less 
computational  burden  in  actual  implementation. 

The  final  model  proposed  at  AFIT  for  use  in  an  extended 
Kalman  filter  was  the  constant  turn-rate  target  acceleration 
model.  Target  trackers  based  on  this  alternative  dynamics 
model  were  initially  developed  for  both  air-to-air  (29)  gun¬ 
nery  (estimation  in  three  dimensions)  and  FLIR  focal  plane 
(5)  target  intensity  tracking  (in  two  dimensions).  Although 
the  extended  Kalman  filter  based  on  a  constant  turn-rate 
target  acceleration  model  was  found  to  outperform  both  an 


extended  Kalman  filter  based  on  a  Brownian  motion  accelera¬ 
tion  model  and  a  multiple  acceleration  model  adaptive  extended 
Kalman  filter  in  two-dimensional  estimation,  subsequent  per¬ 
formance  analyses  (8,21)  indicated  that  comparable  performance 
can  be  achieved  in  the  majority  of  scenarios  through  use  of 
the  enhanced  correlator/linear  Kalman  filter  combination 
based  on  a  first-order  Gauss-Markov  target  acceleration  dynam¬ 
ics  model.  For  cases  where  the  target  being  tracked  is  at 
close  range  or  is  performing  evasive  maneuvers,  the  extended 
Kalman  filter  based  on  the  constant  turn-rate  target  accel¬ 
eration  model  significantly  outperforms  all  of  the  other 
filters  which  have  been  outlined.  However,  the  computational 
burden  imposed  by  a  software  implementation  of  the  constant 
turn-rate  filter  as  previously  defined  (5)  has  strongly 
motivated  use  of  the  linear  Kalman  filter. 

The  computational  burden  is  due  to  two  factors.  First, 
incorporation  of  the  constant  turn-rate  target  acceleration 
model  yields  nonlinear  propagation  equations  as  compared  to 
the  linear  propagation  equations  generated  by  use  of  the 
first-order  Gauss-Markov  target  acceleration  model.  Secondly, 
the  constant  turn-rate  filter  as  originally  defined  used  the 
individual  FLIR  pixel  intensities  as  its  input.  This  design 
required  the  measurement  update  equations  to  be  nonlinear  as 
well,  a  problem  not  encountered  when  the  enhanced  correlator 
produces  target  position  offsets  from  the  center  of  the  FLIR 
field  of  view.  Combining  the  additional  processing  require¬ 
ments  due  to  these  two  sources  of  nonlinearities  with  the 
requirement  of  real-time  performance  compatible  with  the  FLIR 
sensor's  30  Hz  frame  rate,  it  is  easy  to  see  the  preference 
accorded  the  enhanced  correlator/linear  Kalman  filter  ap¬ 
proach. 

However,  the  filter  design  which  yields  the  best  perfor¬ 
mance  against  the  widest  range  of  target  scenarios  is  the 
constant  turn-rate  nonlinear  filter.  Use  of  this  filter  thus 
requires  that  a  method  of  implementation  of  the  algorithm 


having  sufficient  system  throughput  capacity  be  developed, 
that  the  constant  turn-rate  dynamics  model  be  combined  with 
the  linear  measurement  model  of  the  enhanced  correlator,  or 
both.  To  develop  an  optimal  aiming  system  (given  the  devel 
opments  to  date) ,  the  computational  burden  problem  must  be 
overcome . 

1.3  Optical  Processing  Developments 


The  greatest  single  asset  available  through  the  use  of 
optical  architectures  is  the  ability  to  modulate  an  ex¬ 
tremely  high  density  of  data  in  parallel  over  one  or  more 
spatial  dimensions.  Use  of  this  asset  has  tremendously  in¬ 
creased  system  throughput  capacity  over  comparable  electronic 
counterparts . 

An  additional  advantage  of  optical  processing  is  the 
"memorylessness"  of  the  components  used.  Specifically,  this 
advantage  allows  a  large  array  of  functions  to  be  accomplished 
by  a  single  architecture  (e.g.,  matrix-vector ,  matrix-matrix, 
and  matrix-raatrix-raatrix  multiplies  as  well  as  matrix  inver¬ 
sion)  .  Optical  architectures  are  not  locked  into  specific 
functions . 

Recent  work  within  the  field  of  optical  processing  has 
included  the  investigation  and  development  of  solutions  which 
incorporate  algorithms  quite  similar  to  those  outlined  in 
Section  1.2.  Methodologies  for  direct  transfer  of  linear 
algebra  to  optical  architectures  have  also  been  outlined. 

These  results,  when  combined  with  previously  available 
optical  processes,  have  produced  powerful  tools  which  cannot 
be  duplicated  with  a  purely  electronic  configuration. 

Recent  developments  in  transfer  methodologies  have  taken 
a  field  once  limited  solely  to  Fourier-domain  analysis  to 
other  transform  applications  (e.g.,  optical  Mellin  trans¬ 
forms)  ,  subsets  of  which  have  high  applicability  to  problems 
where  the  image  pixel  data  is  statistically  either  station¬ 
ary  or  nonstationary  (9).  Optical  architecture  implementa¬ 
tions  of  such  specific  functions  as  matrix  inversion  (3,4,22), 


matrix-matrix  multiplication  (3.4.10),  matrix-vector  pro¬ 
ducts  (3,4.10),  and  simple  Kalman  filters  (2)  have  been 
detailed  for  bulk  optical  configurations.  Configurations 
for  all  of  these  functions  have  also  been  developed  using 
integrated  optical  architectures  (27),  an  area  which  has 
significant  potential  for  size,  weight,  and  power  consump¬ 
tion  reductions  as  well  as  possessing  an  inherent  advan¬ 
tage  in  system  isolation  stability. 

Furthermore,  use  of  these  architectures  has  been  demon¬ 
strated  for  a  wide  variety  of  applications,  many  of  which  are 
almost  directly  transferable  to  portions  of  the  pointing  and 
tracking  problem.  The  electro- optic  inversion  of  a  complex 
covariance  matrix  generated  by  a  phase-array  radar  (22) 
represents  a  potential  solution,  should  a  phased  array  op¬ 
tical  imaging  configuration  be  ^implemented .  Several  recent 
developments  in  optical  tunable  notch  filtering  (26)  can  be 
reconfigured  for  the  bandpass  instead  of  the  bandreject 
case  and  represent  significant  signal  enhancement  potential. 
Multi-dimensional  correlators  (1)  have  already  been  deployed 
and  are  performing  several  highly  complex  algorithms  in 
real  time.  Although  all  of  these  systems  have  a  large  poten¬ 
tial  for  future  optimization,  their  current  configurations 
as  referenced  represent  an  increase  in  system  bandwidth  over 
their  electronic  counterparts  of  from  two  to  four  orders  of 
magnitude . 

Each  of  these  architectures  can  be  characterized  by  the 
type  of  modulator  used  to  perform  the  processing  operations. 
This  is  due  to  the  fact  that  the  performance  abilities  of  any 
optical  architecture  are  essentially  limited  by  the  operating 
characteristics  (modulation  dimensions,  response  time,  modu¬ 
lation  efficiency,  and  dynamic  range)  of  the  modulator,  since 
it  almost  always  is  the  only  non-passive  system  element  and 
thus  is  subject  to  more  corruptive  processes.  However, 
because  its  role  is  so  critical,  modulator  design  and 


development  have  an  extremly  high  priority  in  the  field  of 
optics,  and  both  many  different  types  and  many  differnt 
modulators  are  available  to  the  optical  architect.  There 
are  several  different  types  of  point  modulators  (not  dis¬ 
cussed  here  because  they  are  generally  not  applicable  to 
optical  processing) ,  as  well  as  linear  and  two-dimensional 
modulators.  Linear  modulators  (30)  are  generally  limited 
to  either  acousto-optic  or  electro-optic  effects  as  the 
type  of  modulation,  whereas  two-dimensional  modulators 
use  many  different  types  of  modulation  criteria.  For  exam¬ 
ple,  included  within  the  subset  of  available  two-dimensional 
modulators  are  liquid  crystal  light  valves  (23),  magneto- optic 
modulators  (4),  and  deformable  mirror  assemblies  (10).  Each 
of  these  types  of  modulators  has  spawned  a  series  of  archi¬ 
tectures  which  are  capable  of  performing  various  functions 
and  which  maximize  the  operational  characteristics  of  the 
particular  device  while  circumventing  its  limitations  as  well 
as  possible.  All  of  the  previously  outlined  architectures 
use  one  or  another  of  these  modulators,  with  each  specific 
problem  solution  having  components  matched  to  the  particular 
modulator  to  produce  an  optimal  solution  insofar  as  the 
optics  is  concerned. 

The  use  of  optical  architectures  will  be  encountered 
throughout  the  research  as  a  method  of  overcoming  the 
limitations  of  the  best  algorithms  currently  available. 


A  more  thorough  review  of  the  optical  developments  used 
will  be  included  within  the  appropriate  chapters. 

1.4  Thesis  Problem 

Laser  weaponry  is  currently  envisioned  for  use  in  a 
multiplicity  of  scenarios:  ground-to-air,  air-to-air,  and 
space- to- space .  As  required  system  performance  as  well  as 
system  operational  parameters  vary  with  application,  limi¬ 
tation  of  the  design  to  one  particular  scenario  is  indi¬ 
cated.  Since  initial  feasibility  studies  and  testing  will 


be  for  a  ground-based  version  of  the  system,  the  problem 
scope  will  be  limited  to  this  scenario.  Additionally,  as  a 
simpler  problem,  it  is  more  useful  for  a  f easibility-of-con- 
cept  demonstration. 

The  proposed  problem,  therefore,  is  the  complete  mathe¬ 
matical  development  and  performance  analysis  of  a  pointing 
and  tracking  system  having  an  optical  data  processor  as  its 
core  for  a  ground-based  laser  weapon.  The  scope  of  the 
problem  will  thus  include  model  development,  choice  of  sen¬ 
sors  exhibiting  optimal  throughput  as  well  as  low  noise 
characteristics,  and  recommendations  for  future  efforts 
based  on  current  developments.  The  mathematical  development 
and  all  performance  analyses  will  be  forcibly  restricted  to 
real  world  constraints  by  limiting  all  component  selection 
to  off-the-shelf  units  having  well-characterized  performance 
criteria. 

The  overriding  goal  of  this  research  will  be  to  demon¬ 
strate  that  the  benefits  offered  by  the  use  of  optics  are 
both  available  and  deployable  now.  Additionally,  the  thesis 
will  demonstrate  that  use  of  optical  techniques  will  enable 
the  employment  of  even  more  sophisticated  dynamics  models 
(should  they  be  developed  in  the  future)  by  having  more  than 
sufficient  system  bandwidth  in  reserve  or  obtainable  through 
alternate  component  selection.  These  demonstrations  will 
form  the  central  core  of  the  proposed  approach. 

1  .5 _ Approach 


The  overall  system  concept  plan  embodies  a  heretofor 
unexplored  tracking  algorithm:  the  enhanced  correlator 
linear  measurement  model  combined  with  the  nonlinear  con¬ 
stant  turn-rate  dynamics  model  as  described  in  Section  1.2. 
While  there  are  several  justifications  for  using  this 
approach  (all  of  which  will  be  outlined),  the  central  fact 
remains  that  the  constant  turn-rate  target  acceleration 
model  exhibits  equal  or  better  performance  over  the  full 


operational  range  of  the  tracker  against  anticipated  targets 
compared  to  any  other  model  yet  developed.  Additionally, 
since  it  has  been  demonstrated  that  use  of  the  constant  turn- 
rate  target  acceleration  model  provides  better  than  adequate 
performance  (as  well  as  equal  of  better  performance  as  noted 
earlier) ,  any  multiple-filter  algorithm  would  merely  compli¬ 
cate  and  add  an  additional  computational  burden  to  the  over¬ 
all  system  without  substantial  performance  benefit.  Such  an 
algorithm  was  considered  for  this  effort  and  consisted  of 
simultaneously  optically  processing  both  the  first-order 
Gauss-Markov  target  acceleration  model  and  the  constant  turn- 
rate  target  acceleration  model  in  parallel.  Both  dynamics 
models  used  the  output  from  an  enhanced  correlator  as  the 
measurement  update.  The  appropriate  dynamics  model  was 
selected  by  considering  the  range  and  range-rate  data  from 
the  previous  and  current  measurement  frames  to  determine  the 
switchover  point.  However,  this  concept  was  rejected  due  to 
the  analysis  outlined. 

The  method  by  which  the  proposed  combination  is  imple¬ 
mented  is  heavily  dependent  on  the  use  of  optics  and  optical 
processors.  Not  only  are  the  Kalman  filter  update  and  propa¬ 
gation  cycle  linear  algebra  operations  accomplished  within 
the  bounds  of  an  optical  processing  architecture,  but  the 
input  signal  is  also  passed  through  a  previously  developed 
(23)  optical  processor  in  order  to  provide  the  cross-corre¬ 
lations  necessary  for  centroid  offset  calculations  at  the 
FLIR  plane.  Use  of  this  technique  thus  allows  for  a  po¬ 
tentially  significant  increase  in  system  processing  band¬ 
width  if  the  FLIR  array  sensor  is  replaced  by  a  real-time 
sensor,  thereby  providing  performance  potential  far  ex¬ 
ceeding  that  obtainable  with  the  FLIR.  This  issue  is  po¬ 
tentially  more  significant  than  any  other  possible  filter 
modifications  and  will  be  discussed  further  in  Chapter  V. 

To  develop  this  admittedly  very  general  overall  descrip¬ 
tion  of  system  operation  into  the  solution  of  the  problem 


As  l  jtr.  Cue  cyr.anies  and  Cue  measurement  update  models 
for  tr.e  proposed  solution  nave  already  been  developed,  the 
initial  task  of  this  thesis  will  be  to  determine  the  designs 
of  the  two  optical  processors,  one  to  perform  the  correlations 
and  one  to  implement  the  filter.  These  designs,  as  noted 
previously,  will  be  constrained  to  tne  real  world  through 
limitations  as  to  component  choices. 

It  is  extremely  important  at  this  point  to  differentiate 
carefully  between  the  two  optical  processors  involved.  This 
is  particularly  true  (as  will  be  seen  in  Chapters  III  and  V) 
since  each  processor  performs  unique  tasks  unrelated  to  the 
other  and  thus  must  be  designed  to  different  specifications. 
For  these  reasons,  the  optical  processor  which  produces  the 
necessary  input  cross-correlations  will  be  termed  the  en¬ 
hanced  optical  correlator  (SOC) ,  whereas  the  optical  pro¬ 
cessor  which  performs  the  filter  operations  will  be  termed 
the  optical  Kalman  filter  (OKF). 

As  the  design  of  the  correlation  processor  has  already 
been  established  (23)  and  will  be  discussed  in  Chapter  III, 
criteria  for  the  design  of  the  OKF  must  be  discussed.  Here, 
choice  of  the  optical  architecture  to  be  used  is  dependent 
on  three  key  parameters:  (1)  the  efficiency  with  which  it 
accomplishes  the  linear  algebra  required  for  Kalman  filter 
operations;  (2)  its  ability  to  perform  the  required  pro¬ 
cessing  without  decreasing  system  dynamic  range  (i.e.,  adding 
noise  to  the  system  as  outlined  further  in  Chapter  V) ;  and 
(3),  its  potential  for  increasing  the  system  throughput  rate. 
Combining  these  requirements  with  the  knowledge  that  any 
optical  architecture  is  limited  by  the  performance  of  the 
modulator,  a  system  tradeoff  analysis  based  upon  the  avail¬ 
able  architectures  can  be  performed. 
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This  tradeoff  analysis  immediately  limits  the  number 
of  architectures  available  for  the  OKF.  This  is  due  to 
several  facts.  First,  since  any  two-dimensional  modulator- 
based  optical  architecture  (designed  for  linear  algebra) 
is  limited  to  a  maximum  system  bandwidth  of  30  Hz  for  pro¬ 
cessing  at  the  full  dynamic  range  of  the  modulator  (4,10,13). 
these  architectures  are  not  suitable  candidates  due  to 
requirement  (3).  Secondly,  electro-optic  modulator-based 
optical  architectures  can  also  be  discounted  due  to  their 
presently  limited  dynamic  range  (30).  This  leaves  only  those 
optical  architectures  based  on  acousto-optic  modulators  as 
suitable  candidates. 

Although  architectures  based  on  acousto- optic  modulators 
can  be  either  one-  or  two-dimensional,  all  of  the  available 
two-dimensional  architectures  require  integrating  detectors 
(1)  and  thus  do  not  satisfy  the  final  requirement  estab¬ 
lished.  This  leaves  only  two  optical  architectures,  both 
based  on  single  acousto-optic  modulators,  as  candidates  for 
the  OKF.  These  are  the  (4)  engagement  array  architecture 
(a  pseudo-parallel  structure)  and  the  (2)  pipelined  itera¬ 
tive  optical  systolic  array  architecture  (a  true  parallel 
structure).  Of  these  two,  the  iterative  nature  of  the  pipe¬ 
lined  iterative  optical  systolic  array  (as  fully  developed 
in  Chapter  V)  proved  far  superior  in  satisfaction  of  the 
efficiency  requirement  as  it  could  accomplish  every  operation 
required  without  major  restructuring.  It  is  for  these  rea¬ 
sons  that  the  pipelined  iterative  optical  systolic  array 
architecture  is  the  solution  of  choice  for  the  OKF. 

Now  that  these  steps  have  been  accomplished,  the  opti¬ 
cal  architecture  for  each  processor  must  be  analyzed  so  as 
to  take  full  advantage  of  the  optical  processing  benefits. 
This  phase  will  thus  require  determination  of  the  required 
system  bandwidths,  choice  of  all  the  specific  optical  com¬ 
ponents  to  be  used,  and  analysis  of  any  constraints  and  the 
methods  by  which  they  are  overcome. 


12 


The  above  steps,  having  been  successfully  cc 
will  then  allow  the  final  design  of  each  optical 
and  computation  of  each  processor's  modulation  tr 
function  (MTF) .  The  MTF,  being  analagous  to  an  i 
response  in  the  electronic  domain,  will  not  only 
system  simulation  and  performance  analyses  from  i 
output,  but  will  also  conclusively  demonstrate  th 
the  optical  architectures  will  not  require  retuni 
filter  due  to  additional  noise. 

If  time  had  permitted,  construction  of  a  lim 
width  breadboard  version  of  each  optical  processo 
have  then  been  accomplished .  Characterisation  of 
response  characteristics  as  compared  to  their  pre 
could  have  been  scaled  upward  and  compared  with  t 
dieted  system  performances  obtained  earlier  for  t 
scale  designs.  Correlation,  identification,  and 
of  observed  deficiencies  in  either  comparison,  if 
would  have  been  reincorporated  into  the  package. 

The  final  step  of  the  research  will  be  to  dr 
sions  pertinent  to  the  proposed  optical  implement- 
the  algorithm  and  to  indicate  areas  of  developmen 
the  potential  to  produce  useful  system  refinement 
that  component  selection  will  be  limited  to  curre: 
able  hardware,  optimizations  in  this  area  may  pro> 
icant  results.  These  results  due  to  hardware  dev 
may  then  impact  the  analysis  within  this  thesis, 
significantly.  Furthermore,  recommendations  as  t 
optical  implementations  of  the  developed  architec 
may  also  hold  vast  potential. 

The  overall  approach  as  outlined  was  selecte< 
provide  as  much  benfit  as  possible  to  both  the  re 
and  the  sponsoring  organization  in  the  five  man-m 
were  available  for  task  completion.  It  is  hoped 
effort  will  be  a  cornerstone  for  future  system  de 


The  definition  of  a  state-space  truth  model  requires 
that  the  mathematical  representation  of  any  real  world  pro¬ 
cess  embodied  within  its  confines  has  the  greatest  degree  of 
correlation  with  the  physical  reality  of  the  real  world  pro¬ 
cess,  given  the  bounds  of  a  finite  state-space.  As  such,  it 
becomes  the  standard  against  which  any  proposed  filter  design 
is  measured  and  thus  mandates  that  the  truth  model  be  free  of 
any  biases  favoring  any  particular  reduced  order  filter  de¬ 
sign  approach. 

In  this  chapter,  a  truth  model  is  developed  for  the 
tracking  portion  of  the  overall  system  for  which  the  design 
goal  is  the  portrayal  of  measured  and  true  motion  of  a  target 
acquired  by  a  ground-based  high-energy  laser  weapons  system. 
The  targets  are  further  postulated  to  be  either  unmanned  mis¬ 
siles  or  manned  aircraft  (single  or  multiple  engine).  To 
perform  this  task  given  the  outlined  constraints,  the  approach 
chosen  (7,8,20,21,23)  has  been  to  project  onto  the  two-dimen¬ 
sional  FLIR  image  plane  the  true  intensity  function  generated 
by  a  target  in  three-dimensional  inertial  space,  translation- 
ally  shifted  to  account  for  the  phase  distortion  caused  by 
the  atmosphere.  Other  physical  phenomena,  such  as  signal 
degradation  due  to  FLIR  optics  mirror  jitter,  which  may  im¬ 
pact  adversely  upon  system  performance  in  a  mathematically 
strict  sense,  are  assumed  to  be  of  negligible  import  com¬ 
pared  to  the  FLIR  image  motion  sources  (due  to  the  true  tar¬ 
get  dynamics  and  the  translational  shift  which  accounts  for 
the  wavefront  phase  distortion  by  the  atmosphere)  and  are 
therefore  omitted  in  the  formulation  of  the  state-space  dy¬ 
namics.  However,  the  simulated  image  that  results  from  the 
described  projection  of  the  target's  true  intensity  function 
onto  the  FLIR  image  plane  is  subjected  to  spatially  corre¬ 
lated  but  temporally  uncorrelated  noise  to  account  for 


background  and  inherent  FLIR  noises  (7).  The  corrupted  image 
simulation  thus  obtained  is  the  measurement  array,  _z(t^),  of 
average  pixel  intensities  (from  an  8x8  pixel  tracking  window) 
that  would  occur  at  the  FLIR  image  plane,  and  is  the  measure¬ 
ment  basis  used  by  the  enhanced  correlator  of  the  filter. 

In  its  original  development,  the  approach  chosen  to  gen¬ 
erate  the  continuous  time  target  model  used  the  stochastic 
form  of  state-space  differential  equations  in  order  to  account 
both  for  the  statistical  properties  of  the  atmospheric  cor¬ 
ruption  and  the  uncertainty  due  to  the  motion  sources  assumed 
to  be  negligible  enough  not  to  warrant  explicit  modeling. 
However,  the  equations  were  developed  in  augmented  form  so 
that  the  states  representing  the  true  target  dynamics  could 
be  mathematically  manipulated  so  as  to  allow  for  deterministic, 
as  well  as  stochastic,  inputs  without  affecting  the  statistical 
properties  of  atmospheric  propagation.  The  approach  thus 
allows  truth  model  generation  of  deterministic  target  trajec¬ 
tories  which,  being  subject  to  all  of  the  corruptive  (and 
statistical)  processes  described  earlier,  can  then  be  used  to 
analyze  the  performance  of  any  proposed  filter  design  approach. 

To  present  this  particular  approach  to  the  design  of  the 
truth  model  in  a  logical  manner,  this  chapter  is  divided  into 
a  series  of  mathematical  descriptions  which  originate  at  the 
true  dynamics  of  a  single  hot  spot  target  and  progress  in  a 
temporal  sense  to  the  resultant  FLIR  image  plane  measurement 
array  used  as  the  input  data  to  the  proposed  filter  design. 
Extensions  for  the  case  of  a  multiple  hot  spot  target  are 
presented  after  complete  development  of  the  single  hot  spot 
target  case.  The  truth  model  thus  developed  can  be  directly 
appended  to  any  proposed  filter  model  design, 

2.2  True  Target  Centroid  Dynamics  Model 

The  projection  of  any  three-dimensional  process  onto  a 
two-dimensional  measurement  plane  requires  an  exact  singular¬ 
ity-free  mathematical  conversion  between  the  two  coordinate 
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systems  m  order  to  be  able  to  describe  the  set  of  all  three- 
dimensional  processes  in  terms  of  the  two-dimensional  coor¬ 
dinate  system.  Since  this  condition  is  fulfilled  within  the 
development  of  this  truth  model,  description  of  the  true  tar¬ 
get  dynamics  centroid  model  is  accomplished  by  mathematically 
specifying  the  four  deterministic  target  test  trajectories 
in  a  three-dimensional  inertial  space  (x^»  y^  ,  and  z^.)  having 
its  origin  at  the  tracker.  This  procedure  allows  specifica¬ 
tion  of  the  conversion  between  the  two  coordinate  spaces  to 
be  developed  at  a  more  logical  point  in  the  overall  truth 
model  development. 

Since  specific  target  trajectories  must  be  developed  not 
only  to  illustrate  the  proposed  filter's  tracking  performance 
under  the  full  range  of  conditions  anticipated  to  be  encoun¬ 
tered  but  also  to  establish  their  baseline  performance,  Tra¬ 
jectory  1  is  defined  as  a  simple  target  traversal  parallel 
to  the  Xj-y^  plane  and  is  depicted  in  Figure  2.  Defining  the 


Figure  2  -  Trajectory  1  Geometry 
inertial  target  position  as 
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Additionally,  the  model  for  Trajectory  1  includes  the  capa¬ 
bility  of  performing  the  traversal  described  with  target  roll 
rates  of  0,  0.5,  and  1  radian  per  second.  This  provision  of 
the  model  allows  inclusion  of  the  change  in  the  FLIR-depicted 
target  shape  as  a  function  of  the  target  profile  as  described 
in  Section  2.5,  and  thus  permits  a  realistic  baseline  to  be 
established . 

Alternatively,  the  definitions  of  Trajectories  2  and  3 
represent  the  desire  to  establish  a  proposed  filter's  track¬ 
ing  performance  against  the  most  extreme  maneuvers  which  any 
potential  target  can  perform.  Trajectory  2,  depicted  in 
Figure  3.  thus  portrays  a  target  which,  after  2.0  seconds  of 
Trajectory  1  motion,  performs  either  a  constant  2g  ( u  =  19.6 

milliradians  per  second)  or  5g  (u>  =  49  milliradians  per 

P 

second)  pull-up  maneuver  over  two  inertial  dimensions  (a 
plane  in  inertial  space)  ending  at  six  seconds  into  the  simu- 


yj (to ) 


X j (to) ,Zj(t) 


Figure  3  -  Trajectory  2  Geometry 

lation.  Since  the  maneuver  start  time  is  at  2.0  seconds 
after  Trajectory  2  initiation,  including  the  results  due  to 
Trajectory  1  motion  for  the  first  two  seconds  yields  the 
mathematical  definitions  of  Trajectory  2  after  initiation  of 
the  maneuver  as 

£j(t)  =  -cos  [cjp  ( t-2)]  km/sec 
y-j(t)  =  sin  [u)p(t-2)]  km/sec 
and  zT(t)  =  0 


(2-6) 

(2-7) 

(2-8) 


Simple  integration  then  yields  the  post  maneuver  target  loca¬ 
tion  as  (based  on  the  target's  initial  inertial  coordinates 
as  defined  by  Trajectory  1) 

Xj(t)  -  5  -  2  -  { sin [u)p( t-2)]  /  km  (2-9) 

y-j(t)  =  0.5  +  ({1  -  cos  [u)p( t-2)]  }  /  tOp)  km  (2-10) 

z];(t)  =  20  km  (2-11) 

An  alternative  form  of  Trajectory  2  is  defined  to  portray 
target  motion  into  the  negative  inertial  z-direction  (resul¬ 
ting  in  the  target  turning  in  toward  the  tracker)  instead  of 
the  inertial  y-direction.  Performance  of  this  trajectory  re¬ 


defines  the  y^  and  equations  above  as 

yx(t)  =  0  (2-12) 

yx ( t )  =  0.5  km  (2-13) 

2j(t)  =  -sin [wp( t-2)]  km/sec  (2-14) 

zT(t)  =  20  -  ({1  -  cos  [u>  ( t-2)]  }  /  u  )  km  (2-15) 

i.  P  P 


with  all  variables  as  defined  earlier.  Use  of  this  modifi¬ 
cation  enables  the  projection  of  three  distinct  ellipsoidal 
target  intensity  patterns  onto  the  FLIR  image  plane  so  that 
target  shape  effects  (to  be  discussed)  can  be  included  in  the 
simulation. 

Noting  that  both  Trajectory  2  descriptions  simulate  the 
trajectory  path  change  as  a  step  function,  it  can  be  conclu¬ 
ded  that  any  proposed  filter  is  required  to  track  a  far 

harsher  environment  than  a  real  situation  where  w  would 

P 

build  up  smoothly  (in  a  continuous  sense)  from  zero  to  its 
specified  terminal  value.  Satisfaction  of  filter  perfor¬ 
mance  specifications  will  therefore  ensure  better  than 
minimal  performance  against  actual  targets. 

The  next  test  trajectory  incorporated  in  the  analysis, 
Trajectory  3,  is  specifically  designed  to  evaluate  the  per¬ 
formance  of  a  proposed  filter  design  against  targets  both 
initiating  and  terminating  a  maneuver.  Implementation  of 
this  trajectory  can  be  considered  as  an  extension  of 


Trajectory  2  in  that  2g  and  5g  puil-up  maneuvers  are  initiated 
and  terminated  as  step  functions  at  2.0  and  3.5  seconds,  re¬ 


spectively.  The  target  has  the  same  initial  position  and 
inertial  velocity  descriptions  as  Trajectory  2  until  maneuver 
termination.  At  t  =  3.5  seconds,  the  target  continues  at  the 
inertial  velocities  (determined  by  use  of  the  Trajectory  2 


equations) 


Xj(t>>3.5)  =  -.99958  km/sec 
y-j.  ( t>3 .5)  =  .029069  km/ sec 
4T(t>3.5)  =  0 


(2-16) 

(2-17) 

(2-13) 


for  the  remainder  of  the  simulation. 


The  final  test  case.  Trajectory  4,  represents  the  por¬ 
trayal  of  target  trajectory  changes  simultaneously  in  all 
three  inertial  dimensions,  rather  than  maintaining  az  least 
one  inertial  dimension  constant  as  in  all  of  the  previous 
trajectory  descriptions.  Use  of  this  test  trajectory  thus 
also  enables  evaluation  of  filter  performance  against  targets 
whose  intensity  pattern  is  continuously  changing  over  all 
three  inertial  dimensions,  which  is  a  more  realistic  por¬ 
trayal  of  the  projected  (assumed)  ellipsoidal  intensity 
patterns  than  that  afforded  by  the  alternate  form  of  Tra¬ 
jectory  2. 

Trajectory  4  is  implemented  by  simulating  a  target 
which,  after  establishment  of  the  baseline  as  defined  for 
Trajectory  2,  performs  a  2g  pull-up  maneuver  which  mathe¬ 
matically  establishes  a  new  three-dimensional  coordinate 
system  as  depicted  in  Figure  4.  For  this  trajectory,  the 


0  =  to  ( t-2 ) 
XI 

T 


Figure  4  -  Trajectory  4  In-  and  Out-of-Plane  Coordinate 

Transformation  Geometry 
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translated  coordinate  system  ( x^, ,  y,p,  and  )  rotates  at  the 
constant  rate  to  about  both  the  xT  and  zT  coordinate  axes. 

The  relationship  between  the  two  coordinate  systems  can  thus 
be  easily  established  by  analyzing  the  velocity  vector  rela¬ 
tionships  for  each  of  the  two  two-dimensional  case  (the  Xj-yT 
plane  and  the  yj~zj  plane)  and  combining  the  results  vec- 
torially.  For  the  trajectory  as  defined,  the  coordinate 
transformations  due  to  the  x^-y^  plane  target  motion  are 

Xj  ( t>2)  =  cos  [to  ( t-2)]  xT(  t)  +  sin  [cop(  t-2)]  yT(  t)  (2-19) 
yI(t>2)  =  -sin  [wp(  t-2)]  xT(  t)  +  cos  t-2)'J  y^(  t )  (2-20) 

Zj ( t>2)  =  ST(t)  *  (2-21) 


where  all  velocity  variables  are  expressed  in  km/sec.  For 
the  out-of-plane  coordinate  case  depicted  in  Figure  4,  the 
coordinate  transformations  due  to  the  y.p-Zrp  plane  target 
motion  are 

xI(t>2)  =  *T(t)  (2-22) 

Vj  ( t>  2)  =  cos[u)p(t-2)]£T(t)  +  sin[ojp(t-2)]  &T(t)  (2-23) 

(t>2)  =  -sin[w  (t-2)]^T(t)  +  cos  [wp(  t-2)]  t)  (2-24) 

where  the  velocity  variable  units  are  defined  similarly  as 
in  the  earlier  case.  Expressing  these  coordinate  transfor¬ 
mations  in  matrix  form  and  combining  them  with  the  initial 
velocity  vector  specified  earlier  (from  Trajectory  2)  yields 
the  overall  coordinate  transformation  vector  after  maneuver 
initiation  as 

Xj(t)  100  cos  0  sin  0  0 

$^(t)  =  0  cos  0  sin  0  -sin  0  cos  3  0 

_2j(t)_  0  -sin  0  cos  8^  __  0  0  1_ 

-cos[ui  (t-2)]  km/sec 

cos  [»  (t-2)]  sin[_Wp(t-2)J  km/sec  (2-25) 

-sin2  t-2)]  km/sec 

where  0  is  defined  in  Figure  4.  Simple  integration  of  the 
overall  coordinate  transformation  then  yields  (in  km) 


(2-26) 

(2-27) 

(2-28) 


xT(t>2)  =  5  -  2  -  sin [jj  (t-2)]/u 
I  p 

y^(t>2)  =  0.5  +  sin2  [uj  ( t-2)  J  /  2u> 

zT(t>2)  =  20  -  (t-2)/2  +  sin[2u  (t-2)]/4w 
X  p 

Finally,  definition  of  the  inertial  coordinate  relationships 
allows  the  out-of- ( inertial) plane  angle  (OPA)  of  the  maneuver 
described  by  the  translating  plane  to  be  easily  generated 
from  the  relationship  (as  seen  in  Figure  4) 

OPA  =  tan"1  aI(t)/yI(t)  (2-29) 

for  all  time  after  maneuver  initiation. 

With  the  true  target  centroid  dynamics  described  in  in¬ 
ertial  space,  the  apparent  FLIR  image  plane  motion  must  be 
transformed  to  the  inertial  motion  descriptions.  Although 
this  transformation  can  be  accomplished  within  several  dif¬ 
ferent  formats,  the  previously  developed  singularity-free 
mathematical  conversion  defines  the  inertial  space  to  FLIR 
image  plane  geometrical  relationship  as  depicted  in  Figures 
5  and  6,  where  the  inertial  coordinates  are  as  defined  ear¬ 
lier,  rg(t)  represents  the  horizontal  tracker  to  target 
range,  and  r(t)  represents  the  tracker  to  target  slant  range. 

Tracker  ,  ^  ■, 

Optics  Port  XI 

(Look  Down  View) 

Zj  ( t) 

2I 

Figure  5  -  Horizontal  Inertial  Space/FLIR  Image  Plane 
Geometrical  Relationship 

It  is  apparent  from  the  figures  that  defining  the  tracker's 
optics  port  as  the  inertial  origin  of  coordinates  results  in 
a  tremendously  simplified  geometry  for  hemispheric  tracking 
while  insuring  the.  singularity-free  conversion  required  for 
the  analysis. 


Target  Centroid 


Figure  6  -  Vertical  Inertial  Space/FLIR  Image  Plane 
Geometrical  Relationship 


The  governing  correlation  equations  are  generated  from 
the  appropriate  figure.  Thus,  in  the  horizontal  direction, 

a(t)  =  tan-"^  [zj  ( t) /x^(  t)]  rad  (2-30) 

yielding  upon  differentiation 

xT(t)aT(t)  -  zT(t)xT(t) 

d ( t)  =  -  rad/sec  (2-31) 

x  j  2 ( t )  +  zI2(t) 

where  Equation  (2-31)  can  be  converted  into  pixels  per  second 
by  using  the  pixel/field  of  view  relationship  of  20  microra¬ 
dians  per  pixel.  Similarly,  in  the  vertical  direction. 

Figure  6  yields  the  equation 

6(t)  =  [y-j- ( t) /r^(  t)]  rad  (2-32) 

where  r^(  t)  =  £xj2(t)  +  z^2(t)J^  as  defined  in  Figure  5. 
Differentiation  thus  yields 


S(t) 


rH(t)^i(t)  -  yx  ( t) f-H( t) 
rH2(t)  +  yI2(t) 
rH(t)yi(t)  -  yj(t)fH(t) 
r2  ( t ) 


rad/ sec 


(2-33) 


where  0(t)  can  also  be  expressed  in  pixels  per  second  through 
use  of  the  same  relationship  as  defined  for  the  horizontal 


direction.  These  equations  coupled  with  the  deterministic 
target  trajectories  (or  a  target  trajectory  expressed  in 
stochastic  form,  as  will  be  discussed  in  Section  2.3)  thus 
fully  define  the  true  target  centroid  dynamics  model  as  pro¬ 
jected  onto  the  FLIR  image  plane.  Furthermore,  definition 
of  the  FLIR  horizontal  and  vertical  axes  in  this  manner  tre¬ 
mendously  simplifies  formulation  of  the  FLIR  image  plane 
target  dynamics  model  which  accounts  both  for  the  target 
motion  and  the  translational  shift  due  to  atmospheric  phase 
distortion.  Therefore,  discussion  of  the  FLIR  image  plane 
target  dynamics  model  is  the  next  logical  link  in  the  model 
chain. 

2.3  FLIR  Image  Plane  Target  Dynamics  Model 

With  both  the  inertial  coordinates  and  the  deterministic 
target  trajectories  defined,  the  model  must  account  for  the 
target  dynamics  as  perceived  by  the  FLIR  image  plane.  The 
approach  chosen  has  been  to  generate  equations  which  repre¬ 
sent  the  motion  due  to  the  target  dynamics  in  terms  of  the 
FLIR  image  plane  coordinate  axes  as  defined  above  (whose 
states  are  subscripted  with  D  for  dynamics)  to  which  are  aug¬ 
mented  atmospheric  states  (subscripted  A)  to  account  for 
the  translational  shift  along  both  axes  due  to  the  atmo¬ 
spheric  jitter  phenomenom. 

The  continuous-time  FLIR  image  plane  target  dynamics 
model  thus  incorporates  two  dynamics  states  to  represent 
the  true  location  of  the  target  centroid  on  the  two-dimen¬ 
sional  FLIR  image  plane.  As  defined  earlier,  these  states 
are  defined  as  a(t),  the  perceived  FLIR  image  azimuth  (hori¬ 
zontal  position),  and  0(t),  the  perceived  FLIR  image  eleva¬ 
tion.  Putting  these  defined  state  variables  in  linear 
differential  equation  form  yields 

xD(t)  =  FD(t)xD(t)  +  BD(t)uD(t)  +  GD(t)wD(t)  (2-34) 

where 
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*D(t)  =  [<S(t)  3  ( t)  ]  1 

d(t)  =  apparent  horizontal  target  velocity  (azimuth 
rate  of  change)  on  the  FLIR  image  plane 
§(t)  =  apparent  vertical  target  velocity  (elevation 
rate  of  change)  on  the  FLIR  image  plane 
Fg(t)  =  the  2x2  state  distribution  matrix 
x^(t)  =  [a(t)  3(t)]^  with  Qt  ( t )  and  0(t)  as  defined 

above 

Bp(t)  =  the  2x2  control  input  distribution  matrix 
Up(t)  =  the  two-dimensional  vector  of  system  control 
inputs 

Gj-j(t)  =  the  2x2  dynamics  driving  noise  distribution 
matrix 

and  Wjj(t)  =  the  two-dimensional  dynamics  driving  noise 

vector  assumed  to  be  independent  white  noise 
of  mean  zero  and  strength  ^(t);  i.e., 

E{ wD(t) Wq( t+x) } =^D( t) 6 (t) 

The  nonlinear  form  of  the  continuous- time  FLIR  image  plane 
target  dynamics  model  is  defined  to  be 

*D(t)  =  f [xD(b) ,uD(t) ,t]  +  GD(t)wD(t)  (2-35) 

where 

f  ( t)  »Jip  ( t)  ,  t]  =  the  two-dimensional  vector  composed 

of  possibly  nonlinear  functions  of 
Xj-j ( t )  ,  uD(t.)  ,  and  t  which  represents 
the  state  description  of  target 
motion  driven  by  the  inputs  over 
time 

is  the  only  other  definition  required.  Definition  of  the 
stochastic  state-space  differential  equations  in  this  manner 
therefore  allows  the  generation  of  the  true  target  trajecto¬ 
ries  that  are  dictated  by  test  requirements  and  outlined  in 
Section  2.2.  This  deterministic  form  is  accomplished  by  using 
the  linear  form  Equation  (2-34)  while  defining 


.  - .  • .  ■  .  ■  .  -  •  .  ■ 


ZD(tJ  =  0 
BD(t)  =  I 
and  Gp(t)  =  _0 

thereby  allowing  truth  model  generation  of  specific  time 
histories  of  a(t)  and  3 ( t )  as  a  function  only  of  the  deter¬ 
ministic  inputs,  but  in  the  stochastic  differential  equation 
format  by  developing  the  appropriate  a(t)  and  0(t)  to  yield 
these  time  histories  via  integration.  It  is  important  to 
note,  however,  that  redefinition  of  F^(t)  and  G„(t)  to  allow 
nonzero  entries  permits  incorporation  of  either  a  mixed 
deterministic/stochastic  form  (for  only  one  FLIR  axis 
affected)  or  a  fully  stochastic  form  of  Equation  (2-34)  into 
the  truth  model  trajectory  description.  This  capability  can 
be  used,  for  example,  to  simulate  jinking  or  turbulence 
effects  on  top  of  a  purposeful  (deterministic)  trajectory. 

Equations  (2-34)  and  (2-35)  are  then  augmented  by 
expressions  which  represent  the  apparent  target  motion  and 
statistics  which  are  due  to  the  atmospheric  wavefront  phase 
distortion.  These  atmospheric  disturbances,  as  developed 
previously  (7),  are  modeled  as  a  third-order  Gauss-Markov 
process  which  is  describable  as  the  output  of  a  shaping 
filter  with  a  frequency  domain  transfer  function  of 


aA  K( 1 4. 1 4) (659. 5)2 

W.J  ( s+ 1  4 . 1 4)  (  s+ 65q  =  5 ) 2 


(2-36) 


driven  by  a  unit  strength  white  gaussian  noise  w^.  A  dupli¬ 
cate  independent  model  is  also  used  to  generate  3A  in  the 
vertical  FLIR  image  plane  direction.  The  form  of  the  sto¬ 
chastic  differential  equation  used  to  represent  the  atmo¬ 
spheric  jitter  thus  becomes 

jtA(t)  =  FA(t)xA(t)  +  GA(t)wA(t)  (2-37) 


where 
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state 


x^(tj  =  vector  of  the  six  atmospheric  noise 
derivatives 

F^(t)  =  atmospheric  plant  matrix 
G^(t)  =  atmospheric  noise  uncertainty  distribution  matrix 
and  w^(t)  =  two-dimensional  vector  of  white  Gaussian  noise 
inputs  with  statistics:  E{w^(t)}  =  _0  and 
E{wA(t)wJ^t+x)}  =^A(t)6(x) 

After  augmenting  the  atmospheric  states  to  the  dynamics  state 
the  discretized  solution  to  the  truth  model  propagation  (7) 
for  the  motion  of  the  target  centroid  as  perceived  by  the 
FLIR  image  plane  (using  a  deterministic  model  for  true  cen¬ 
troid  motion)  becomes 


"B  ,  ( t .  )]  f  0 

(t^,)  =  i(t4i1  ,t.  )x(t,  )  +  -a-i-  ud(t.)  + -  (2-38) 


i+V  i+1  '  i'-v  “i' 


where 


x(  ti' 


state  vector  of  the  two  dynamics  states  and  the  six 
atmospheric  states 


— D  ^i  +  1  '  ^i^ 

1  0 

0 

1  *A(t1+1.ti). 

exp  FqA t  | 

0 

.  o  T 

exp  F.AtJ  At  = 

-  t . 
i 


B  ,(t.  ) 
— d  i 


-d(tF 


since  both  F^  and  F^  are  time  invariant  for  this 
problem 

input  matrix  for  dynamics 
i  ±( t . . 1 , t . )3( x)dx 

piecewise  constant  function  (between  sample  times) 
evaluated  at  the  midpoint  of  the  interval  as  an 
approximation  to  the  integral  of  d(t)  and  S(t) 
from  t^  to  t^^^ 

<S(t^  +  At/2) 

5 ( t .  +  At/2) 


—Ad 


(t4;  5 


zero  mean  discrete- time  white  Gaussia 

T 


statistics : 

a 


E(«Ad(ti,aAd(tj)) 


15  . 


Ad^f  '  ^  SAdltiJ 


7  iiv 


=  ti.Miti+VT'l*k(x'l*Vzl±A 


T 

!'t)g!  ( 


S{,9-£ 


'i  +  1 


Ad(tiJ  ^Ad(ti))  =  S 


&Ad(ti^ -Ad^i^-Ad^j  ^  ^Ad(tj^  T 
where  ^  is  the  Cholesky  square  r: 

sition  of  (15) 

and  the  apparent  target  centroid  is  located  at  ( 

Explicitly,  the  state  transit: 


[3D(ti)+3A(ti)]). 


[$(ti+^,t^)  from  Equation  (2-38)]  is 
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where  the  lower  two  3x3  submatrices  in  £(t^+^,t^) 
and  propagate  the  atmospheric  disturbances  in  the 
and  elevation  directions  with  the  single  pole  (re 
Pi  above)  and  double  pole  (p2  above)  as  defined  j 
(2-36).  Finally  for  the  purposes  of  this  truth  n 
is  defined  to  be 


Definition  01  all  factors  as  identified  above  thus  allows 
truth  model  representation  of  the  target  dynamics  on  the  FLIP, 
image  plane  with  all  major  factors  affecting  target  motion  in¬ 
cluded.  However,  as  the  actual  FLIP  image  plane  presentation 
consists  of  an  IR  intensity  function  distribution  measurement 
rather  than  a  focused  image  of  the  target  centroid,  the  truth 
model  must  account  for  these  additional  factors  inherent  to 
the  real  world  process.  Definition  of  these  factors  in  the 
FLIR  image  plane  measurement  model  presented  below  thus  com¬ 
pletes  the  truth  model  for  the  single  hot  spot  case,  from 
which  extensions  can  be  made  to  define  the  truth  model  com¬ 
pletely  for  the  additional  case  of  a  multiple  hot  spot  target. 

2.4  FLIR  Ima?e  Plane  Measurement  Model 


In  order  to  preserve  the  maximum  degree  of  correlation 
between  the  truth  model  and  the  physical  reality  of  the  actual 
process  within  the  bounds  of  a  finite  state-space,  any  mathe¬ 
matical  distribution  used  to  represent  the  IR  image  distri¬ 
bution  on  the  FLIR  image  plane  must  approximate  the  actual 
measured  pattern  to  a  very  high  degree.  Previous  efforts 
(7,16,20)  in  this  area  have  shown  that  a  more  than  adequate 
relationship  exists  when  the  pattern  due  to  distant  targets 
is  represented  by  a  bivariate  Gaussian  function  having  cir¬ 
cular  equal  intensity  contours.  Furthermore,  it  has  been 
additionally  established  (7,20)  that  the  only  change  required 
for  representation  of  closer  range  targets  is  the  geometry 
used  to  represent  the  equal  intensity  contours.  For  this 
case,  an  elliptical  pattern  maintains  the  desired  correla¬ 
tion.  These  developments  allow  the  perceived  FLIR  image 
plane  intensity  function  for  either  case  to  be  expressed  as 

l(a,0)  =  IMax  exp{-^[(a-ap)  (0-0p)]  [p_1]  [(a-a?)  (0-0p)]T)  (2-39) 

where  the  variables,  as  shown  in  Figure  7,  are 

IM  =  maximum  target  intensity 


coordinates  of  I.,  which  correspond  to  the 

.'lax 

apparent  target  centroid  location  on  the  FLIP, 
tracking  window 

eigenvalues  of  P  corresponding  to  the  ellipse 
semimajor  axis  along  the  velocity  vector  and 
the  ellipse  semi mi nor  axis  perpendicular  to 
the  velocity  vector,  respectively,  where  1  and 
s  denote  large  and  small  (measured  in  pixels) 
target  velocity  component  perpendicular  to  the 
line  of  sight  from  the  FLIP,  image  plane  to  the 
target,  which  is  defined  to  be  conincident 
with  the  ellipse  semimajor  axis  as  depicted 
v  -|-FLIR  plane  orientation  angle  defined  by 
the  intensity  pattern  semimajor  axis  and  the 
FLIR  horizontal  axis  (a)  translated  to  the 
apparent  target  centroid  (a  ,8  ) 


Figure  7  -  Image  Intensity  Characteristics 


Prior  to  any  additional  development,  one  additional  fac¬ 
tor  must  be  considered.  Recalling  that  each  pixel  of  the  8x3 
FLIR  tracking  array  has  a  field  of  view  (FOV)  of  2 0x20  micro¬ 
radians,  it  is  apparent  that  the  small  resultant  total  FOV  is 
conducive  to  approximating  the  target  angular  displacement  by 
the  linear  displacement  of  the  target  image  centroid  on  the 
FLIR  plane.  It  has  also  been  demonstrated  (7)  that  this  re¬ 
lationship  holds  for  the  FLIR  image  plane  velocity  components 
as  well.  Since,  as  noted  previously,  all  of  the  displacement 
and  velocity  FLIR  plane  components  (a,  4,  6,  and  3)  are 
expressible  in  terms  of  angular  or  displacement  units,  an 
initial  angular  calculation  is  required  to  generate  the  angu¬ 
lar  displacement  of  the  inertial  target  velocity  (depicted  in 
Figure  8)  out  of  the  FLIR  image  plane.  The  entire  balance  of 
the  truth  model  target  intensity  pattern  generation  can  then 
be  accomplished  using  angular  measurements. 

The  first  step  in  this  generation  is  thus  to  calculate 
the  target  velocity  plane  to  FLIR  image  plane  angle,  y.  This 
computation  is  necessary  to  obtain  the  FLIR  image  plane  velo¬ 
city  component  perpendicular  to  the  line  of  sight  (v  ^) ,  which 
is  the  projection  of  the  inertial  velocity  (v)  known  from  the 
simulation.  From  Equations  (2-31)  and  (2-33),  4  and  5,  re¬ 
spectively,  are  generated  in  terms  of  radians  per  second. 

Since  r(t)  is  also  defined  by  Equation  (2-33),  it  is  easy  to 
see  than  an  equivalent  representation  of  FLIR  image  plane 
velocity  translated  to  inertial  space  is 

4'(t)  =  r(t)4(t)  meters/second  (2-40a) 

and  S’(t)  =  r(t)3(t)  meter s/ second  (2-40b) 

where  the  prime  superscript  is  used  to  indicate  variables 
translated  from  the  FLIR  plane  to  the  target  and  whose  dis¬ 
placement  variables  are  always  expressed  in  meters  instead  of 
pixels.  Since  the  inertial  velocity  is  known  to  be 


and  the  FLIP,  image  plane  velocity  (as  depicted  in  Figure  8) 


is  defined  to  be 


vpl  |  =  [<s!(t)  +  i52(t)]J 


(2-42) 


the  target  velocity  plane  to  FLIP  image  plane  angle  y  is 


cos  y  =  r(t) | vpll /|v| 


(2-43) 


The  remainder  of  the  angular  relationships  depicted  in  F'.gures 
8  and  9  can  now  be  developed.  Specifically,  these  are 


cos  0  =  <S(t)/|  v  ^ 
sin  9  =  k{  t)  /  I  v 


(a _,8_)  \&rZ2) _ 


^  pi 


(2-44) 

(2-45) 


Figure  8  -  Image  Projection  Relationships 

With  cos  y  determined,  the  semimajor  axis  of  a  single 
hot  spot  target  can  be  projected  onto  the  FLIR  image  plane. 
For  5  defined  to  be  the  length  of  the  target  in  meters,  the 
apparent  length  of  the  target  to  be  measured  by  the  FLIR 
plane  is 

1  =6  cos  y  meters  (2- 

as  depicted  in  Figure  9.  Since,  as  noted  previously,  the 
fields  of  view  are  sufficiently  small,  the  angle  'V  is 

f  =  aT'/r(t)  radians  (2- 


Figure  9  -  Ellipse  Semimajor  Axis  Projection  Relationships 

which  then  allows  to  be  calculated  as  (with  the  appro¬ 
priate  pixel  conversion  factor) 

=  l{//20x10~k  pixels  (2-48) 

Similarly,  the  ellipse  semiminor  axis  can  be  defined  as 

as  =  Og 1 /r(20x1 0"^)  pixels  (2-49) 

where  1  is  the  target's  IR  cross-section. 

In  this  analysis,  an  additional  manipulation  is  incor¬ 
porated  to  lessen  computation  time.  Recalling  that  each  tra¬ 
jectory  as  defined  yields  y=0  at  t=0,  o^  and  o g  immediately 
become 

0^(0)  =  5/r (20x1 O-^)  pixels  (2-50a) 

Og(0)  =  Og 1 /r(20x1 0~^)  pixels  (2-50b) 

To  preclude  recalculation  of  all  quantities  at  each  analysis 
point,  the  simple  ratios 

aL(t)  =  (cos  y)  [oL(0)]  L'r(0)]/r(t)  ,  t^O  (2-51a) 

as(t)  =  [as(0)]  [r(0)]/r(t)  ,  t^O  (2-51b) 

are  set  up  to  calculate  the  ellipse  semimajor  and  semiminor 
axis  lengths  efficiently.  This  then  establishes  all  param¬ 
eters  required  for  the  analysis. 

The  intensity  at  any  point  in  the  image  plane  can  now  be 
computed  in  the  image  ellipse  principal  axis  coordinate  sys¬ 
tem.  Here,  Equation  (2-39)  is  redefined  after  the  coordinate 


change  as 


0 


-1 


IU.S)  =  IMax  exp{  -  k  [AL  AS] 


AL 

AS 


(2-52) 


where 


AL  =  (a-a  )cos  '0  +  (3-3  )sin  6 
P  P 

AS  =  (3-3  )cos  6  -  (a-a  )sin  8 
P  P 

and  0  is  as  defined  earlier.  Since  the  average  intensity  for 
any  pixel,  as  measured  by  the  FLIR,  is  the  integral  of  the 
apparent  target  intensity  function  over  the  pixel  area,  divi¬ 
ded  by  the  area  of  the  pixel,  and  corrupted  by  FLIR  and  back¬ 
ground  noise,  these  effects  must  be  included  in  the  final  ex¬ 
pression. 

To  generate  this  final  expression,  it  can  be  noted  that 
previous  work  (7)  has  shown  that  the  apparent  target  intensit 
function  can  be  well  approximated  by  averaging  the  intensity 
over  25  equally  spaced  points  within  each  pixel.  Additional 
work  done  at  the  same  time  documented  the  existence  of  spa¬ 
tial  correlations  of  background  noise  in  each  data  frame  with 
nonnegligible  spatial  correlations  between  each  pixel  and  its 
closest  two  neighboring  pixels  in  each  direction.  Since  this 
work  also  justified  that  any  temporal  correlations  of  these 
noises  are  negligible,  only  a  single  noise  term  (n^)  need  be 
incorporated  into  the  expression  for  the  average  intensity  in 
the  pixel  in  the  1-th  row  and  m-th  column  of  the  8x8  array. 
This  expression  now  becomes 


=  I/Max  exPf-iOhjk  aSDjlJ 
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(2-53) 


and  represents  the  measurement  array  for  a  single  hot  spot 
target  which  is  then  input  to  the  correlator/filter.  For 
multiple  hot  spot  targets,  three  Gaussian  hot  spots  with 
elliptical  constant  intensity  contours  and  parallel 


semimajor  axes  were  used  in  this  study.  The  measurement  thus 
becomes  the  summation  of  the  apparent  intensity  functions  due 
to  each  spot  plus  the  noise  term  for  each  pixel.  This  devel¬ 
opment  is  fully  outlined  in  the  following  section. 

2 .5  Pro.iectio n  of  M u  1 1 iple  Hot  Spots 


For  previous  truth  model  calculations  involving  only 
single  hot  spot  targets,  the  target  center  of  mass  (COM) 
was  assumed  conincident  with  the  center  of  the  intensity 
ellipsoid.  Since  the  method  for  projecting  the  target  in¬ 
tensity  pattern  is  the  same  for  both  cases,  a  solution  for 
the  target  COM  for  the  multiple  hot  spot  case  can  potentially 
be  generated  given  the  assumption  that  the  angular  orientation 
of  any  of  the  hot  spots-'  principal  axes  is  identical  to  the 
others.  If  the  relationship  of  each  ellipsoid  center  to  the 
target  COM  is  known  in  inertial  coordinates  as  well,  a  sim¬ 
ple  transformation  of  these  relationships  into  FLIR  image 
plane  coordinates  will  guarantee  a  solution. 

Prior  to  defining  these  additional  relationships  for  the 
thesis  problem,  it  is  necessary  to  set  up  the  geometry  pro¬ 
perly  for  their  use.  As  depicted  in  Figure  10,  where  all 


e0,emC 


e  ,  e 
a  ma 


Figure  10  -  Relation  of  the  eraa~em3  Plane  Unit  Vectors  to  the 
FLIR  Plane  Unit  Vectors 

the  coordinate  frames  are  shown,  the  unit  vectors  e  and  e  , 

ma  m3 

are  related  to  the  FLIR  plane  unit  vectors  (for  ease  in 
implementing  the  required  coordinate  transformation)  bv  the 


expressions 
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^ma 


e 


roi3 


(2-54a) 
(  2  -  5  4  b ) 


where  the  subscript  m  is  used  to  define  FLIP,  plane  translated 
variables.  Figure  11  thus  depicts  the  translation  of  the 
plane  determined  by  e  „  and  e  _  in  inertial  soace  so  that  its 
origin  is  coincident  with  the  target's  COM.  Since  this  new 
plane  is  parallel  to  the  FLIP  plane,  any  descriptions  of  tar¬ 
get  properties  in  terms  of  this  "m"  plane  are  easily  relat- 
able  to  FLIR  plane  variables  as  well. 


From  the  trajectory  model,  the  velocity  vector  of  the 
target  is  known.  Since  the  location  of  hot  spots  in  a  fixed 
frame  can  be  determined  for  a  specific  simulated  target,  the 
velocity  vector  (which  is  assumed  to  point  towards  the  nose 
of  the  target)  can  be  used  to  define  one  axis  of  another 
coordinate  frame  used  to  accomplish  this  hot  spot  determina¬ 
tion  also  having  its  origin  at  the  target  COM.  Defining  the 
unit  vector  in  this  axis  as  e  ,  two  other  unit  vectors  are 
defined  as  follows.  By  performing  a  cross-product  of  e  with 
a  unit  vector  (e.)  in  the  inertial  y-direction,  a  second  axis 
with  unit  vector  e^  which  will  always  lie  in  the  horizontal 
inertial  plane  is  defined.  The  final  axis  of  the  new  frame 
is  then  defined  as  having  a  unit  vector,  e^,  corresponding  to 


v  v'v  '  •  -V/vV' 
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the  cross-product  of  the  two  previously  defined  target  frame 

vectors.  These  relationships  are  depicted  in  Figure  12. 

e  (along  velocity  vector) 
emS  t? 


\ 


Figure  12  -  e  -e  0  and  e  -e  -e  Relationships 
&  ma  mg  x  y  z 

A  coordinate  system  for  the  target  must  now  be  defined. 
Assuming  that  the  target  is  a  multi-engine  aircraft  with  hot 
spots  as  depicted  in  Figure  13.  the  arbitrary  initial  specif 
cation  for  this  thesis  is  that  the  aircraft's  initial  orien¬ 
tation  (of  the  plane  defined  by  the  aircraft  wings  and  fuse¬ 
lage)  is  coincident  with  the  e^-e^  plane.  Additionally,  the 
following  unit  vectors  used  to  establish  a  coordinate  refer¬ 
ence  frame  for  the  aircraft  are  defined:  e  from  the  air- 

-  ax 

craft  COM  towards  the  nose,  e  „  from  the  aircraft  COM  toward 

the  right  wing,  and  e  as  the  cross-product  of  e  and  e 

&  az  ^  ^  ^  ax  ay 


Figure  13  -  Initial  Ellipsoid  Centers 

These  assumptions  and  definitions  thus  specify  that  the  cen¬ 
ter  or  the  ellipsoid  due  to  hot  spot  1  in  Figure  13  will 
always  lie  along  the  e  axis.  Furthermore,  the  centers  of 


the  ellipsoid  due  to  hot  spots  2  and  3  will  remain  along  the 
e _  axis  unless  the  aircraft  performs  a  roll  maneuver. 

For  this  analysis,  roil  maneuvers  were  simulated  by 
applying  one  of  the  constant  roll-rates  as  defined  for  Tra¬ 
jectory  1  (Section  2.2)  as  a  step  input.  In  a  coordinate 
frame  sense,  the  necessary  transformations  are  depicted  in 
Figure  14*  Although  a  step  input  application  of  the  maneuver 
is  again  done  for  simulation  simplicity,  the  technique  in¬ 
sures  that  any  performance  analysis  exhibiting  adequate  or 
better  performance  can  reasonably  be  expected  to  reflect  a 
design  which  will  perform  even  better  in  a  real  world  case 
where  the  roll-rate  would  build  up  smoothly. 


e  ,  e 
x  ax 


Roll  Rate 


y  Tay 

Figure  14  -  Roll  Maneuver  Geometry 

Mathematically ,  the  only  calculation  resulting  from  per¬ 
formance  of  a  roll  maneuver  is  a  simple  planar  rotation  to 
determine  the  direction  of  the  ellipsoid  centers.  For  the 

aircraft  roll-rate  defined  to  be  w  ,  the  amount  of  planar 

r 

rotation  [with  the  assumption  that  t}>(0)=0  and  t  defined  as 
the  time  of  roll  maneuver  initiation]  becomes 


<t>  ( t)  =  u  >  ( t  -  t  ) 
r  r 


( 2-55. 


The  planar  rotation  depicted  in  Figure  14  can  thus  be  ex¬ 
pressed  as 


e  =  e  cos  q  +  e  sin  <t> 
ay  y  *  z  ■ 


2-  56 ) 


and  the  location  of  the  three  hot  spots  relative  to  the  co¬ 
ordinate  frame  specified  by  the  aircraft  velocity  vector  and 


tne  horizontal  inertial  plane  can  always  be  determined. 

All  that  remains  to  be  done  to  complete  the  analysis  is 
to  convert  the  rotated  distance  to  the  equivalent  offset 


distance  on  the  FLIR  image  plane.  If  A 


ay 


.s  defined  to  be 


the  distance  in  meters  from  the  aircraft  COM  to  the  center 


of  hot  spot  2  (lying  in  the  positive  e  ,  direction),  dot 

a  y 

products  of  this  vector  length  with  the  "m"  plane  unit 


vectors  yield  the  offsets  with  the  required  geometry. 

plane  yields 

axes  as 


Specifically,  this  projection  to  the  eraa~em|3 


A  and  A  0 
ma  m3 
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( 2- 57a) 
( 2- 57b) 


as  shown  in  Figures  15  and  16. 
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Figure  15  -  Hot  Spot  Offset 


To  convert  A  and  A  a  to  FLIP,  image  clane  offsets  re- 
ma  m3  ° 

quires  only  simple  calculations  at  this  point.  Since  the 

-+■  -f 

emct_emS  P^ane  as  parallel  to  the  FLIP  plane  and  normal  to  th 
line- of - sight  vector  to  the  aircraft,  the  small  angle  approx 
imation  is  again  valid,  yielding  the  FLIR  plane  displacement 

A  =  0.  =  A  /r(20x10  pixels  (2-58a 

ct  ma  / 

and  A0  =  H7  =  A  o/r(20x10  )  pixels  (2-5Sb 

p  mp 

These  displacements  are  then  at  the  FLIR  plane  location  of 

the  center  of  the  ellipsoid  due  to  hot  spot  2.  Since  hot 

spot  3  is  located  on  the  same  (aircraft)  axis  in  a  reverse 

orientation  and  is  assumed  to  be  at  the  same  distance  from 

the  aircraft  COM  as  hot  spot  2,  its  offset  can  readily  be 

calculated  to  be  ( -A  ,  -A0). 

a  3 

To  demonstrate  this  projection,  previous  efforts  (7,3,2 
have  run  a  simple  trajectory  having  the  initial  inertial  co¬ 
ordinates 

XI  ( 0)  =  yj(0)  =  z-j-(O)  =  10  km 

and  the  constant  inertial  velocities 

Xj(t)  =  y-j(t)  "  ^(t)  =  -0.5  km/sec 

Since  this  defines  a  target  moving  directly  towards  the  cen¬ 
ter  of  the  FLIR  image  plane  and  having  a  velocity  vector  nor 
mal  to  the  FLIR  image  plane, 

& ( t)  =  3 ( t)  =  0 

for  the  entire  trajectory.  Hot  spots  were  defined  as 

1  at  one  meter  forward  of  the  aircraft  COM 

along  e 
&  ax 

2  and  3  at  ±0.5  meters  on  the  sides  of  the  aircraft 
COM  (along  ©ay) 

and  two  simulations  were  used  to  demonstrate  the  projection 
model . 

The  first  simulation,  depicted  in  Figure  17,  shows  that 
for  a  roll  of  tt  radians,  the  ellipsoid  center  due  to  hot  spo 
1  stays  centered  in  the  FLIR  plane,  while  the  ellipsoid 
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centers  due  to  hot  spots  2  and  3  roil  through  tt  r 
while  remaining  at  a  relative  position  (to  each  o 
radians.  The  pixel  offset  change  in  the  figure  ( 
seen  in  Figure  18  which  depicts  a  roil  of  2tt  radi. 
to  the  increase  in  projected  hot  spot  offset  dist. 
function  of  the  decreasing  range.  For  example,  f 
simulation  time  of  five  seconds,  the  initial  pixe. 
for  the  two  hot  spots  in  the  horizontal  FLIR  dime: 


L3x(10,000)2J 5(20x10-d) 
whereas  at  the  termination  of  the  simulation,  the, 


±0.5 


L'3x(  7500)  2]  5(20x10 


1  .93  pi: 


For  a  target  performing  a  roll  through  2tt  rac 
gure  18  clearly  illustrates  the  above  relationshij 
maintaining  the  correct  offsets  as  defined.  Thus 
is  valid  and  can  be  used  for  the  multiple  hot  spol 
Now  that  the  truth  model  used  to  input  data  1 
ter  model  has  been  defined,  the  next  step  of  this 
quires  that  the  tracking  algorithm  to  be  used  be  1 
scribed  prior  to  depicting  the  results  of  the  sysl 
formance  analyses.  Since  a  hybrid  ( optical/elect: 
processor  has  been  determined  to  be  the  solution  c 
for  this  effort,  the  most  logical  method  of  acconn 
this  step  is  to  follow  the  flow  of  data  from  inpul 
As  the  input  stage  requires  a  complete  definition 
enhanced  correlator  and  the  (optical)  method  of  il 
mentation ,  Chapter  III  is  devoted  solely  to  this  1 


III.  Enhanced  Correlator 


3.1  Overview 

The  sole  thrust  of  this  chapter  is  to  define  the  en¬ 
hanced  optical  correlator  (EOC)  completely.  To  accomplish 
this  task  requires  that  the  original  electronic  implementa¬ 
tion  (23)  of  the  correlator  be  both  mathematically  specified 
and  justified  as  well  as  the  optical  processor  to  be  used. 

The  one  additional  analysis  required  to  complete  this  chapter 
is  the  comparison  of  the  projected  EOC  performance  with  that 
obtainable  from  a  software  implementation  after  all  EOC  com¬ 
ponents  have  been  specified. 

Two  points  require  clarification  before  proceeding  with 
the  analysis.  First,  the  software  processed  version  of  this 
algorithm  uses  the  raw  intensity  data  obtained  from  each  FLIR 
frame  as  the  input  to  the  correlator.  The  correlator  then 
processes  this  data  to  obtain  the  noise  corrupted  target 
position  offset  measurements  from  the  center  of  the  FLIR 
tracking  array  field  of  view,  _z  ( t  ^ )  .  In  contrast,  the  opti¬ 
cally  processed  enhanced  correlator  us^s  the  projected  IR 
image  (that  would  be  seen  by  a  FLIR  array  sensor  if  one  were 
present)  as  its  input.  It  then  projects  the  required  corre¬ 
lations  for  a  simple  target  center  of  mass  calculation  onto 
a  measurement  array  sensor.  As  a  FLIR  array  sensor  was  man¬ 
dated  for  these  efforts,  such  a  .se.nsor  is  used  in  all  analy¬ 
ses.  Chapter  V  suggests  some  alternatives  in  this  area 
which  have  the  potential  for  significantly  enhancing  the 
performance  of  the  overall  tracking  system. 

Secondly,  successful  demonstration  of  the  optically  pro¬ 
cessed  enhanced  correlator  thus  requires  that  it  output 
equivalent  or  better  results  than  those  obtainable  from  a 
software  implementation  of  the  algorithm.  This  is  particu¬ 
larly  true  due  to  the  method  by  which  each  version  generates 
the  required  Fourier  transforms  outlined  in  Section  3.2. 
Although  previous  analyses  (23)  have  shown  that  the  algorithm 


-7 


used  in  the  software  version  is  essentially  insensitive  to 
wordlength,  it  is  obvious  that  the  quality  of  the  transforms 
(particularly  those  due  to  images  corrupted  by  noise)  is 
dependent  on  the  specific  discrete  Fourier  transform  (DFT) 
software  used  in  the  correlator.  When  the  requirement  of 
processing  within  the  30  Hz  FLIR  frame  rate  (especially  when 
combined  with  the  additional  processing  required)  is  combined 
with  the  infeasibility  of  having  a  large  main-frame  computer 
at  every  weapon  site,  it  can  be  postulated  that  the  quality 
(in  terms  of  correlation  peak  locations)  may  suffer  in  an 
actual  implementation  due  to  use  of  a  more  efficient  but 
less  accurate  DFT  algorithm.  In  an  optical  version,  the 
limitation  in  transform  production  (at  least  for  the  opti¬ 
cally  produced  transform  of  the  noise-corrupted  signal)  is 
dus  to  two  factors.  Since  analog  processing  is  used,  the 
transform  quality  is  initially  limited  by  the  dynamic  range 
of  the  modulators  used  (here,  approximately  40  dB  at  best). 
Additionally,  since  the  optical  transform  is  produced  spa¬ 
tially  and  optical  components  do  not  have  infinite  dimensions 
some  transform  components  are  lost  (this  will  be  addressed 
in  Sectiqn  3.3).  However,  due  to  the  nature  of  the  algorithm 
outlined  in  Section  3.2,  direct  comparison  of  the  optically 
processed  correlator  to  the  software  processed  version  as 
previously  defined  (23)  is  achievable,  because  the  limita¬ 
tions  outlined  for  the  optical  package  are  of  negligible 
significance  for  data  of  very  low  frequencies  as  in  this  case 
Thus,  the  analysis  in  Section  3.4  provides  a  direct  compar¬ 
ison  of  the  EOC  implementation  to  a  large  main-frame  (Cyber 
175)  implementation. 

With  these  points  clarified,  the  mathematics  to  be  used 
must  be  specified.  The  following  section  performs  this  task 
for  the  algorithm  itself.  Those  interested  in  details  of  the 
software  implementation  may  wish  to  consult  Appendix  A  of 
this  thesis  which  is  identical  to  the  same  section  of  Refer¬ 
ence  23. 
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3.2  Enhanced  Correlator  Algorithm  Development 


The  concept  behind  the  enhanced  correlator  algorithm 
represents  a  desire  to  combine  the  benefits  offered  by  the 
relatively  low  computational  burden  of  a  correlation  algo¬ 
rithm  (as  compared  to  a  highly-dimensioned  extended  Kalman 
filter)  with  the  statistical  accuracy  (in  an  environment 
which  is  not  free  of  sources  of  error)  available  from  a 
Kalman  filter.  Additionally,  the  developed  algorithm,  by 
exploiting  knowledge  about  the  target's  dynamics  and  any 
disturbances  which  have  the  potential  to  cause  apparent 
translational  intensity  function  offsets,  enhances  perfor¬ 
mance  over  standard  correlation  algorithms  which  use  only 
the  raw  data  from  previous  frames. 

To  accomplish  these  actions,  the  template  to  be  used 
for  correlation  is  generated  from  data  obtained  directly 
from  the  propagation  cycle  of  the  Kalman  filter  which  will 
be  described  in  Chapter  IV.  Specifically,  the  propagated 
state  estimates,  x(tT+^),  are  used  to  generate  [from  Equa¬ 
tion  (2-39)J  an  estimated  FLIR  target  intensity  function, 
h  [x(tT  +  .| )  »t^+^]  .  This  estimated  FLIR  target  intensity  func¬ 
tion  has  two  important  properties:  (1)  it  is  due  to  a  target 
center  of  mass  (COM)  located  exactly  at  the  center  of  the 
FLIR  image  plane  since  the  propagation  cycle  of  the  Kalman 
filter  produces  the  necessary  data  to  center  the  target  in 
the  tracking  window  field  of  view  (FOV);  and  (2),  it  is 
uncorrupted  by  either  temporal  or  spatial  noise  by  design. 
These  properties,  combined  with  the  template's  incorporation 
of  image  intensity  changes  due  to  target  roll  rates  and/or 
range  rate  changes  (see  Section  2.5),  would  thus  correlate 
perfectly  with  an  intensity  pattern  due  to  a  perfectly 
modeled  target  in  a  noise-free  environment. 

Since  it  is  as  impossible  to  model  the  target  parameters 
of  interest  perfectly  in  the  real  world  as  it  is  to  have  a 
noise-free  environment,  correlation  of  a  template  represen¬ 
ting  an  assumed  perfect  target  intensity  in  an  ideal  (free 


of  noise)  environment  with  the  signal  due  to  the  actual  pro¬ 
cess  in  a  real  environment  will  yield  the  best  correlation 
that  can  be  expected.  In  either  a  purely  electronic  or  opti¬ 
cal  configuration,  this  correlation  is  best  transformed  in 
the  transform  domain.  Defining 

F[g(a,S)]  =  G(fa,fg)  (3-1) 

and  F[l(a,3)]  =  L(fa.fg)  (3-2) 

to  be  the  Fourier  transforms  of  the  uwo-dimensional  pixel 
sequences  _g(ot,3)  and  _l(a,8),  their  cross-correlation  R(a,3) 
can  be  written  as 

R(a,3)  =£(oi,3)  *l(o,3)  =  F" 1  [G ( f , f g)  •  L*(fa,f3>3  (3~3) 

where  L*(fQ,fg)  represents  the  complex  conjugate  of  L(fQ,fg). 
If  it  is  further  postulated  that  l(a,3)  represents  the  esti¬ 
mated  Fill  target  intensity *f unction  h £x( 1 7+ ^ , t^+ ^  )J  from  the 
Kalman  filter  and  _g(ot,3)  represents  the  actual  FLIR  target 
intensity  sequences,  it  is  easy  to  see  that  implementation  of 
the  algorithm  of  Appendix  A  yields  the  required  Fourier 
transforms  required  for  the  enhanced  correlator. 

To  examine  the  correlation  in  specific  detail,  consider 
the  simple  case  of  correlating  two  two-dimensional  Gaussian 
structures.  Here,  l(a,3)>  as  depicted  in  Figure  19,  is  the 
template  due  to  a  single  hot  spot  target  which  is  sufficiently 
far  enough  away  from  the  tracker  (or  a  target  with  a  circular 
cross-section  having  a  velocity  vector  normal  to  the  FLIR 
plane)  to  have  circular  equal  intensity  contours  with  a2  =  2 
pixels.  Furthermore,  the  representation  of  the  target  as 
depicted  reveals  that  the  template  data  is  padded  by  sur¬ 
rounding  the  original  8x8  array  by  eight  rows  and  columns  of 
zeros  to  generate  a  24x24  pixel  array  structure.  This  arti¬ 
ficial  bandwidth  limitation  (as  the  Fourier  frequencies  are 
now  present  in  a  spatial  sense)  of  an  assumed  periodic  in 
two  dimensions  mathematical  function  which  in  reality  is  in¬ 
finite  in  extent  allows  the  FFT  algorithm  to  represent  the 
transform  far  more  accurately  than  for  the  case  of  a  non -pad¬ 
ded  array  where  higher  frequencies  are  undefined. 
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Figure  19  -  Centered  Single  Gaussian  Template 

Similarly,  Figure  20  depicts  a  typical  representation  of 
£(ot,8).  This  array  of  slightly  noise-corrupted  data,  con¬ 
taining  the  target's  offset  intensity  function,  is  offset  by 
one  pixel  in  the  horizontal  (eQ)  direction  from  the  center  of 
the  24x24  pixel  array.  Correlation  of  the  data  depicted  in 
Figures  19  and  20  then  produces  the  cross-correlation,  R(a,3)» 
depicted  in  Figure  21 . 

Noting  that  the  expected  result  of  correlating  two 
Gaussian  arrays  is  another  Gaussian  array,  it  can  easily  be 
established  that  the  symmetry  of  the.  correlation  array  will 
be  representative  of  any  offsets  (spatial  frequency  differ¬ 
ences)  present  between  the  two  input  data  arrays.  The  data 
in  Figure  21  shows  this  information  explicity.  In  the  verti¬ 
cal  (e^)  direction,  R(a,@)  is  symmetric  which  indicates  that 
no  offset  is  present  in  this  dimension  between  the  target  and 
template.  The  fact,  however,  that  a  one  pixel  circular  shift 
to  the  left  is  required  to  restore  symmetry  (a  function  of 
the  assumed  origins  of  the  two  images  in  performing  the  DFT) 
in  the  horizontal  direction  clearly  shows  that  R(a,g)  does 
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Figure  20  -  Noise  Corrupted  Data  Array 
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Figure  21  -  Result  of  Cross-Correlation  of  Data  and  Template 


indeed  contain  the  _g(a,3)  offset  information  due  to  the  tar¬ 
get’s  true  position  relative  to  the  tracking  window  at  the 
sample  time. 

Before  proceeding  further  with  the  analysis,  it  is  im¬ 
portant  to  note  that  use  of  the  circular  shift  is  justified 
by  the  assumption  that  the  Fourier  transforms  produced  are 
in  reality  infinitely  periodic.  Since  the  data  contained  in 
R(cj,3)  represents  only  a  single  period  of  the  result,  the 
next  pixel  to  the  left  of  column  one  is  defined  to  be  equal 
to  the  corresponding  pixel  in  column  24.  Similarly,  corre¬ 
sponding  pixels  in  rows  one  and  24  are  defined.  Once  these 
definitions  are  made,  the  correlation  data  depicted  in  Fi¬ 
gure  21  can  be  restructured  to  a  more  useful  form. 

This  R(a,0)  data  restructuring  is  depicted  in  Figure  22 
Here,  in  line  with  the  definitions  established  above,  the 
original  correlation  data  (Figure  21)  of  quadrant  I  is 
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Figure  22  -  Diagonal  Quadrant  Swap  of  Cross-Correlation 


swapped  with  the  original  correlation  data  of  quadrant  IV. 

In  a  similar  manner,  the  data  of  quadrants  II  and  III  are 
also  swapped.  This  180-degree  circular  shift  over  two  dimen 
sions  produces  a  R(a,0)  representation  where  the  relative 
position  offset  between  the  template  and  target  is  now  the 
offset  of  the  resulting  Gaussian’s  maximum  from  the  center 
of  the  sequence.  This  form  thus  enables  an  easy  and  effi¬ 
cient  calculation  of  the  offsets  as  will  be  described. 

It  is  important  to  note  at  this  point  that  this  shift 
is  only  required  for  a  purely  electronic  implementation  of 
the  algorithm.  This  is  due  to  the  fact  that  the  origin  of 
the  correlation  depicted  in  Figure  21  directly  corresponds 
to  the  origins  of  the  original  two  data  sequences  depicted 
in  Figures  19  and  20.  Use  of  this  convention  requires  that 
the  origins  of  the  original  two  sequences  be  superimposed  to 
create  the  first  correlation  data  point.  Thus  the  correla¬ 
tion  data  corresponding  to  the  four  pixels  which  surround 
the  centers  of  the  FOVs  of  the  original  two  sequences  is 
contained  in  the  four  corner  pixels  of  the  correlation 
matrix.  By  defining  the  center  of  the  24x24  correlation 
matrix  as  the  correlation  of  the  centers  of  the  original 
sequences,  the  enhanced  optical  correlator  directly  pro¬ 
duces  the  data  depicted  in  Figure  22. 

Returning  to  the  analysis,  a  comparison  of  the  original 
Gaussian  template,  _l(a,0),  depicted  in  Figure  19  with  the 
R(a,8)  output  shown  in  Figure  22  reveals  some  degeneracies 
which  must  be  addressed.  Specifically,  these  degeneracies 
are  a  spreading  of  the  function  as  well  as  some  degree  of 
assymmetry  caused  by  the  noise  present  in  the  _g(a,0)  se¬ 
quence.  Since  the  R(a,0)  pixel  magnitudes  correspond  di¬ 
rectly  to  the  degree  of  resemblance  between  the  noise-free 
template  and  the  noisy  target  signal,  simple  thresholding 
of  the  R(a,0)  data  can  produce  significant  smoothing  of  the 
correlation.  The  use  of  thresholding  is  justified  if  the 
assumption  is  made  that  the  correlations  produced  by  the 
noise  and  template  are  insignificant  in  magnitude  compared 


to  the  correlations  of  interest.  Additionally,  the  use  of 
thresholding  all-  s  an  easy  calculation  of  the  COM  rather 
than  a  more  complicated  and  ambiguous  peak  detection  scheme. 
Figure  23  thus  depicts  the  result  of  setting  the  threshold  a 
half  the  value  of  the  maximum  magnitude  of  the  correlation 
shown  in  Figure  22.  The  data  of  Figure  23  was  produced  by 
subtracting  the  threshold  value  from  each  pixel  in  turn.  If 
the  resulting  magnitude  was  less  than  or  equal  to  zero,  the 
element  was  considered  to  contain  poor  correlation  informa¬ 
tion  and  thus  was  set  to  zero.  All  resulting  positive  raagni 
tudes  were  retained  at  their  original  ( pre- thresholding) 
values.  As  can  clearly  be  seen  in  Figure  23,  significant 
smoothing  of  both  degeneracies  is  induced. 
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Figure  23  -  Result  of  Thresholding  of  .5 

The  fraction  of  the  sequence  maximum  magnitude  used  as 
the  thresnold  warrants  careful  investigation.  If  it  is  set 
too  low,  the  correlations  due  solely  to  the  noise  are  not 
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significantly  countered  and  the  resulting  ouM  statistic 
suffer.  The  original  development  (23)  of  the  enhanced  cor¬ 
relator  used  0.5  as  the  threshold  fraction.  Further  efforts 
in  this  area  (21)  have  indicated  a  potential  for  enhanced 
performance  if  a  value  of  0.3  is  used.  Thus,  although  the 
example  used  for  illustrative  purposes  above  uses  the  0.5 
threshold  value,  the  analyses  of  the  enhanced  correlator 
operation  as  defined  in  this  thesis  directly  represent  the 
use  of  the  lower  stated  value. 

Now  that  the  data  has  been  obtained  in  a  usable  form, 
an  efficient  means  of  extracting  it  is  required.  Instead  of 
a  peak  detector,  the  enhanced  correlator  uses  a  centroid 
summation  of  the  R(a,$)  sequence  since  it  is  assumed  that 
the  center  of  mass  of  the  threshold  correlation  is  a  good 
indication  of  the  peak  location.  Specifically,  for  either 
the  horizontal  or  the  vertical  direction,  this  centroid 
summation  is  defined  as 


where  i  is  the  horizontal  or  vertical  coordinate  of  a  given 
pixel,  ll^l  represents  the  magnitude  of  the  pixel  intensity 
at  that  location,  and  N  is  the  total  number  of  pixels  in  the 
array.  Equation  (3-4)  is  thus  used  twice  in  the  correlator, 
one  time  to  produce  the  horizontal  centroid  and  the  second 
time  to  produce  the  vertical  centroid.  The  resulting  cen¬ 
troid,  jC(a,3),  is  then  the  measured  position  offset  of  the 
target  from  the  center  of  the  data  frame. 

The  development  o-f  Equation  (3-4)  concludes  the  speci¬ 
fication  of  the  enhanced  correlator.  Since  this  thesis  en¬ 
visions  the  use  of  an  optical  architecture  to  perform  this 
enhanced  correlation,  the  design  of  the  enhanced  optical 
correlator  must  be  specified  and  its  performance  analyzed. 
The  next  section  thus  specifies  the  SQC  design. 


Prior  to  EOG  design  specification,  an  outline  of  the 
basic  laws  of  Fourier  optics  as  they  pertain  to  this  develop¬ 
ment  is  required.  Since  this  outline  can  hardly  be  termed 
comprehensive,  the  reader  interested  in  exploring  this  topic 
further  is  urged  to  consult  Reference  6. 

Starting  with  the  simple  case,  consider  first  the  system 
illustrated  in  Figure  24.  This  two-dimensional  system  uses 
(as  indeed  do  all  optical  transform  systems)  optical  inter¬ 
ference  to  process  incoming  signals.  Specifically,  light 
from  a  laser  point  source  S  is  collimated  by  lens  L  .  This 
collimated  beam  of  quasi-monochromatic  coherent  light  (as 
defined  further  in  Chapter  V)  then  illuminates  an  object 
placed  in  plane  Pi  resulting  in  the  formation  of  a  diffrac¬ 
tion  pattern.  By  placing  the  object  exactly  one  focal  length 
( f i )  in  front  of  the  transform  lens  Li,  a  phase-flat  classi¬ 
cal  Fraunhofer  diffraction  pattern,  t0(x2,y2),  due  to  the 
object's  space- varying  amplitude  transmittance  is  formed  at 
plane  P2  located  at  exactly  one  focal  length  behind  Li.  The 
complex  field  of  this  Fraunhofer  diffraction  pattern  is 
mathematically  equal  to  the  spatial  Fourier  transform  of  the 
object.  This  can  clearly  be  illustrated  by  comparing  both 
mathematical  definitions,  which  are 


To(x2,y2)  =  kj / /t o ( xi , y i ) exp 


I  Afi 


■(x2xi  +  y2y  i ) 


dxidy i 


To  ( f  v  >  f =  ki//to  (xj  ,yi)  expf- j2m(f  xi  +  f  yi )  ]  dxidyi 

a  y  x  j 

where 


(3-5) 

(3-6) 


ki  s  a  complex  constant  required  for  each  integral 
To(x2,y2)  =  Fraunhofer  diffraction  pattern  at  P2 
To(fx»f  )  =  Fourier  transform  of  the  spatial  frequencies 
of  the  object  at  Pi  at  P2 
and  A  =  wavelength  of  the  laser 

As  the  quantities  (x2/Afi)  and  (y2/Afi)  do  indeed  define  th 
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horizontal  and  vertical  spatial  frequencies,  resj 
due  to  the  object  at  Pi,  the  two  expressions  are 
tically  identical. 

With  a  Fourier  transform  available  at  ?2,  me 
of  its  information  can  be  simply  produced.  By  p] 
optical  filter  at  P2,  the  complex  optical  filter 
produces  a  second  diffraction  which  is  then  coll? 
L2.  Defining  f2  to  be  the  focal  length  of  L2,  a 
image  is  produced  at  plane  P3  when  configured  as 
Figure  24. 

Extending  these  principles  to  the  dynamic  cs 
EOC,  simple  extensions  from  the  static  case  outii 
serve  to  define  the  required  process  completely, 
image  formed  at  the  FLIR  image  plane  is  a  dynamic 
of  the  static  object  specified  at  Pi.  Since  this 
actually  _g(oc,3),  it  can  easily  be  seen  that  G(f  , 
produced  at  P2.  Since  this  transform  must  be  mul 
L(fa,fg)  in  order  to  produce  the  required  correla 
R(ot,3),  some  type  of  modulator  is  required  to  ”fi 
G(fa,fg)  data  points  on  a  real-time  (or  pseudo  re 
basis.  It  is  important  to  note  that  this  modulat 
able  to  represent  L(fa,f^)  as  defined  in  Section 
suming  for  the  moment  that  such  a  modulator  is  av 
it  is  obvious  that  the  resulting  diffraction  patt 
reality  the  multiplication  of  the  two  transforms 
by  Equation  (3-3) •  Thus,  as  L2  is  in  effect  prod 
Fourier  transform  of  a  Fourier  transform,  the  res 
verted  image  at  P3  is  R(a,3).  By  placing  an  arra 
at  this  image  plane,  the  required  cross- correlati 
readily  be  converted  back  into  an  electronic  sign 
electronic  signal  is  then  used  to  compute  _C(a,3) 
thresholding  all  of  the  data. 

Since  the  type  of  array  detector  to  be  used 
been  specified,  the  only  true  component  variable 
specified  is  the  modulator  to  be  used.  Although 


the  requirement  of  proper  lens  specification,  the  statement 
is  valid  due  to  the  fact  that  modulator  choice  will  specify 
lens  choice,  since  the  modulator  dimensions  will  determine 
the  spacing  between  the  discreet  transform  frequencies  which 
a  lens  must  provide.  Here,  lens  choice  is  completely  deter¬ 
mined  by  the  relative  spatial  dimensions  of  the  FLIR  array 
sensor  pixels  and  the  modulator  pixels  as  the  FLIR  image 
plane  (Pi)  image  is  assumed  to  correspond  exactly  with  the 
FLIR  array  sensor.  Thus,  an  examination  of  the  performance 
requirements  of  the  modulator  will  serve  to  specify  the 
appropriate  choice  of  modulator  from  those  available  and 
the  balance  of  the  system  as  well. 

Looking  at  the  obvious  requirements,  the  modulator  must 
both  be  able  to  modulate  the  IR  signals  of  interest  effi¬ 
ciently  and  be  able  to  perform  this  modulation  at  a  minimum 
system  throughput  rate  of  30  Hz.  The  efficiency  requirement 
arises  from  two  factors:  (1)  the  IR  signals  of  interest 
cover  a  fairly  wide  dynamic  range  (at  least  30  dB  given  the 
original  use  of  the  FLIR  sensor  at  the  image  plane) ;  and 
(2),  not  all  modulators  are  designed  to  operate  at  the  IR 
wavelengths  of  interest.  Additionally,  detectors  at  these 
wavelengths  (30)  do  not  exhibit  high  quantum  efficiencies 
(the  ratio  of  electrons  emitted  per  photons  incident). 

Since  the  minimum  system  throughput  rate  precludes  the  normal 
method  of  longer  integration  times  to  produce  larger  signals, 
the  efficiency  requirement  is  the  determining  factor.  Ad¬ 
ditionally,  as  all  of  the  modulation  devices  outlined  in 
Chapter  I  are  capable  of  operation  at  30  Hz,  the  efficiency 
requirement  becomes  the  mandatory  performance  criterion 
factor  as  well. 

Other  performance  criteria  exist  which  would  enhance 
(or  not  detract  from)  system  performance  and  thus  are  desir¬ 
able.  Foremost  among  these  is  pixel  dimension  compatibility 
(with  the  FLIR  array  sensor  pixel  size)  so  that  extensive 
image  magnification  which  can  lead  to  optical  aberrations  can 
be  avoided.  Secondly,  an  ability  to  support  a  system  throughput 
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rate  greater  than  30  Hz  (while  maintaining  the  modulation 
efficiency)  is  desired.  Finally,  since  the  Fourier  transform 
L(fa,fg)  is  produced  in  software,  direct  digital  loading  to 
the  device  is  desirable  to  avoid  the  delay  incurred  by  D/A 
conversion.  Basically,  these  desired  performance  criteria 
indicate  a  desire  to  provide  a  noise-free  correlator  (des¬ 
pite  thresholding  elimination  of  any  errors  that  would  be 
produced)  in  terms  of  the  optics  themselves  as  well  as  pro¬ 
viding  the  ability  to  support  higher  bandwidth  tracker  opera¬ 
tions  . 

Out  of  all  the  modulators  outlined  in  Chapter  I,  all  of 
these  criteria  merge  in  only  one  distinct  modulator.  This 
modulator,  manufactured  by  Texas  Instruments,  is  the  TI 
deformable  mirror  assembly.  From  the  TI  specification  sheet 
(Number  82-149,  March  1982),  the  performance  criteria  of 
interest  are:  % 

Modulation  Bandwidth  120  Hz 

Dynamic  Range  >42  dB  from  DC  to  100  Hz 

Modulation  Efficiency  at 

Wavelengths  of  Interest  0.88 

Input  Signal  Required  14-bit  Digital 

Pixel  Dimension  Same  as  FLIR  array 

Pixel- to-Pixel  Accuracy  ±1  bit 

Since  all  specifications,  both  mandatory  and  desired,  are  met 
or  exceeded  only  with  this  modulator,  it  is  obvious  that  this 
is  the  modulator  of  choice.  It  must  be  noted,  however,  that 
choice  of  this  modulator  requires  thermoelectric  coolers  for 
optimal  performance.  This  is  due  to  the  fact  that  the  mirrors 
on  the  array  are  metallic  and  will  partially  absorb  the  IR 
signals.  Since  this  requirement  is  easily  met,  the  final  de¬ 
sign  can  now  be  specified. 

As  illustrated  in  Figure  25,  this  design  is  quite  simple. 
It  is  assumed  both  here  and  in  the  analysis  of  Section  3.4 
that  optical  bandpass  filtering  to  produce  an  image  of  only 
the  signals  of  interest  is  accomplished  by  the  input  optics. 
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Figure  25  -  Enhanced  Optical  Correlator  Design 

The  choice  of  lenses  illustrated  is  a  direct  result  of  po¬ 
sitioning  the  optical  transform  products  G(fa»fg)  produced 
by  Li  at  the  centers  of  each  pixel  of  the  deformable  mirror 
array.  Additionally,  both  Li  and  L2  have  the  same  focal 
length  so  as  to  maintain  the  same  image  size  at  both  input 
and  output  planes.  The  only  other  pertinent  fact  of  inter¬ 
est  is  the  option  made  available  by  the  larger-tnan-required 
dynamic  range.  This  option  allows  for  image  intensification 
in  the  following  manner.  Since  the  magnitude  of  I,j  is 
known  from  the  previous  frame  of  data,  the  next  frame  can 
boost  or  lower  the  maximum  value  (and  thus  all  other  values) 
of  L(fa,fg)  to  compensate  at  least  partially  for  low  dynamic 
range  signals.  Use  of  this  option,  although  unexplored  for 
this  effort,  would  thus  maintain  FLIR  sensor  operation  at  as 
efficient  a  level  as  is  possible. 

Now  that  the  EOC  design  has  been  specified,  an  analysis 
of  its  performance,  particularly  compared  to  an  enhanced 
correlator  implemented  totally  in  software,  is  required  for 
system  acceptance.  This  process  is  accomplished  next  in 
Section  3.4* 


Before  developing  the  system  MTF,  the  system  operating 
characteristics  must  be  specified.  For  the  purpose  of  this 
thesis,  these  characteristics  are  specified  to  be:  (1)  sys¬ 
tem  operation  at  a  throughput  rate  of  30  Hz;  (2)  input  of 
the  maximum  intensity  of  L(fa,fg)  at  +30  dB  to  the  modulator; 
and  (3),  use  of  the  FLIR  array  sensor  as  defined. 

The  EOC  MTF  is  developed  here  only  to  illustrate  the 
flatness  (i.e.,  the  uniform  response)  of  the  optical  system. 
Since  the  image  cannot  be  altered,  the  only  parameters  of 
interest  are  the  system  components.  With  the  lens  diameters 
specified  to  be  at  least  four  times  the  image  size  at  the 
image  plane,  the  lenses  produce  no  significant  aberrations 
that  would  affect  image  quality  as  defined.  Similarly,  use 
of  a  FLIR  array  sensor  does  not  produce  any  noise  which  has 
not  been  previously  modeled.  This  leaves  only  the  defor¬ 
mable  mirror  array  as  a  potential  source  of  error. 

Since  this  modulator  is  specified  as  having  ±1  bit 
accuracy,  the  worst  case  is  considered.  Assuming  two  con¬ 
tributing  (at  the  image  plane  P3)  one-bit  errors  are  made 
while  operating  at  the  30  dB  (10-bit)  dynamic  range  speci¬ 
fication,  the  relative  error  would  be 

%  variation  =  (2/1024)  x  100  =  0.195*  (3-7) 

which  is  insignificant.  Thus  the  image  produced  at  the  FLIR 
array  sensor  is  subject  to  a  maximum  error  that,  if  consid¬ 
ered  noise,  would  be  far  below  the  other  noise  sources  of  the 
problem . 

The  major  analysis  of  interest  is  a  comparison  of  the 
optical  transform  quality  with  that  obtainable  from  a  soft¬ 
ware  implementation.  Here  the  parameter  of  interest  is  the 
relative  spatial  frequencies  obtainable  by  each  method.  If 
the  optical  system  can  produce  transforms  of  equal  or  better 
resolution  than  those  implemented  in  software,  use  of  the 
EOC  is  enabled.  For  this  case,  the  limiting  system  resolu¬ 
tion  will  be  the  transform  developed  in  software,  L(f  ,f,). 

—  a  p 


Since  this  transform  is  limited  by  the  FLIR  array  sensor 
(or  modulator)  spatial  resolution  in  either  the  EOC  or  a 
software  implementation,  the  center-to-center  pixel  distance 
becomes  the  limit.  As  this  distance  is  specified  (by  the 
same  specification  sheet  outlined  earlier  in  this  section) 
to  be  one  millimeter  for  either  dimension,  the  maximum  spa¬ 
tial  resolution  required  is  one  per  millimeter.  As  the  field 
of  Fourier  optics  can  easily  generate  spatial  resolutions  of 
100  per  millimeter  (6),  the  spatial  resolution  required  is  well 
within  the  achievable  bounds  of  the  EOC  and  its  use  is  jus¬ 
tified  . 

With  R(a,0)  available,  software  will  be  used  to  generate 
C(a,(3)  after  thresholding.  This  process  is  described  in  de¬ 
tail  in  Chapter  V.  Since  C ( ot ,  0 )  is  the  desired  measurement 
to  be  incorporated  within  the  filter,  the  next  logical  step 
in  the  development  of  the  proposed  tracker  is  to  specify  the 
filter  models  used.  Only  after  the  filter  algorithm  has  been 
described  can  the  optical  Kalman  filter  be  specified  and  the 
entire  system  subjected  to  performance  analyses.  Thus,  the 
next  chapter  is  devoted  to  this  task. 


IV.  Filter  Models 


Inherent  to  any  Kalman  filter  algorithm  is  the  strict 
separation  of  its  two  distinct  functions,  propagation  and 
measurement.  Since  the  Kalman  filter  is  designed  to  operate 
in  a  stochastic  (uncertain)  environment  as  well,  each  filter 
design  requires  two  models  given  the  assumption  that  the 
measurement  uncertainty  is  truly  independent  of  the  uncer¬ 
tainty  due  to -incomplete  knowledge  of  the  system's  dynamics. 

The  first  of  these  required  models  is  the  measurement  update 
model,  which  is  used  to  generate  enhanced  state  (variable  of 
interest)  estimates  from  measurements  which  are  assumed  to 
be  corrupted  by  noise.  The  other  required  model  propagates 
the  enhanced  state  estimates  forward  in  time  based  on  know-  % 
ledge  of  the  system's  dynamics  (which  are  assumed  tp  be 
known  with  some  degree  of  uncertainty).  As  such,  it  is 
termed  the  system  dynamics  model. 

Although  a  relatively  new  field  in  filtering,  a  large 
number  of  stochastic  filter  algorithms  are  available  (11,15,17). 
These  include  both  linear  and  nonlinear  Kalman  filters  as 
well  as  several  types  of  second-order  filters,  some  of  which 
can  be  implemented  either  on  a  continuous- time  or  a  sampled- 
data  basis.  Here,  the  sampled-data  basis  is  required  by  the 
problem  definition.  However,  the  problem  definition  does  not 
specify  a  choice  as  to  the  type  of  stochastic  filter  algorithm 
to  be  employed. 

Noting  that  the  measurement  update  and  sytem  dynamics 
models  for  this  problem  have  been  implemented  in  both  linear 
and  nonlinear  form  in  the  past,  a  set  of  tradeoff  criteria 
needs  to  be  developed  to  make  a  realistic  design  choice 
(given  the  necessity  of  operating  within  the  FLIR  30  Hz  frame 
rate)  from  the  available  options.  As  in  all  other  engineering 
professions,  the  most  important  criterion  is  the  ability  of 
the  system  to  perform  at  an  acceptable  level  over  the  full 
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range  of  the  problem  specification.  However,  particularly  in 
military  systems,  better  than  the  defined  acceptable  perfor¬ 
mance  level  is  often  desired.  This  may  necessitate  accepting 
the  additional  computational  burdens  of  a  nonlinear  Kalman  or 
second-order  stochastic  filter  algorithm  in  order  to  provide 
more  realistic  modeling.  If  such  is  the  case,  the  system 
sampling  rate  may  have  to  be  decreased  due  to  the  additional 
computation  time  incurred.  Since  fewer  measurements  are  thus 
available  for  a  given  time  period,  the  performance  enhance¬ 
ment  obtained  may  not  be  justifiable.  This  is  particularly 
true  if  the  decrease  in  sampling  rate  requires  a  major  re¬ 
structuring  of  previously  specified  system  comnnents. 

As  outlined  in  Chapter  I,  this  thesis  problem  arises 
from  this  dilemma.  Here,  two  methods  of  solution  are  used 
to  generate  an  achievable  filter  algorithm  (in  terms  of  the 
time  specification)  with  enhanced  performance:  (1)  optical 
processing  techniques  with  their  inherently  higher  processing 
rates  are  used  extensively  throughout  the  system;  and  (2), 
a  linear  measurement  model,  deemed  to  have  sufficient  accu¬ 
racy  for  the  desired  tracking  system  performance,  is  used  in 
place  of  a  nonlinear  algorithm  in  an  attempt  to  reduce  the 
overall  computational  burden.  By  using  both  of  these  solu¬ 
tions,  the  disadvantages  of  additional  computational  burden 
are  overcome  while  yielding  the  performance  enhancement  due 
to  more  realistic  modeling  of  the  problem's  dynamics. 

To  define  the  filter  in  the  most  logical  manner,  sepa¬ 
rate  sections  are  devoted  to  each  filter  model.  In  a  tempo¬ 
ral  sense,  the  next  system  activity  after  production  of  the 
centroid  coordinates,  C(a,0),  is  the  Kalman  filter  measure¬ 
ment  update  cycle.  As  such,  the  next  section  fully  describes 
this  model.  Section  4.3  then  concludes  this  chapter  with  a 
description  of  the  filter's  propagation  cycle  as  specified 
by  the  state  dynamics  model. 


4.2  Measurement  Update  Model 


Recalling  that  the  input  _C(a,3)  to  the  measurement  up¬ 
date  cycle  of  the  filter  is  in  reality  the  target's  center- 
of-mass  (COM)  position  in  the  tracking  window,  a  linear  law 
is  easily  developed.  By  defining 

C  ( a ,  3)  =  [aQ  3c]T  (4-1) 

where  and  3c  are  the  measured  target's  COM  pixel  positions 
at  a  given  sample  time  t^ 

2(t.)  =  [oc(ti)  3c(ti)]T  (4-2) 

where  _z(t.)  is  the  discrete- time  measurement  vector  required 
by  the  model. 

Before  specifying  the  model,  however,  several  other 
definitions  are  required.  First,  the  state  vector  itself 
must  be  specified.  Here  an  eight-state  state  vector  is 
defined  as 


where  the  subscript 
superscript  denotes 


an 


A  A 

0D  %  v6 
denotes  a 
estimate 


AAA 

ao  a 
ot  3 

filter 


A 


variable , 


(4-3) 

the  caret 


^D’  E  ^*IR  plane  target  position  estimates  due  to 
the  target's  dynamics 

vQ ,  Vg  =  FLIR  plane  target  velocity  estimates 

a^ ,  a.g  =  FLIR  plane  target  acceleration  estimates 

and  ,  3^  =  FLIR  plane  target  position  estimates  due  to 

atmospheric  effects  on  the  signal  of  interest 


Additional  variables  requiring  definition  are 

tT  =  time  immediately  prior  to  measurement  avail¬ 
ability,  to  which  all  state  and  state  uncer¬ 
tainty  values  are  propagated 
tt  =  time  postulated  to  be  immediately  post  mea¬ 
surement  which  represents  state  estimate  and 
state  covariance  (uncertainty)  measurement-en¬ 
hanced  values 
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v(t,.  )  =  additive  noise  corruption  at  the  sample  time, 
assumed  to  be  white,  and  Gaussian,  with  sta¬ 
tistics  to  be  determined 

R(t^)  =  the  degree  of  measurement  uncertainty  due  to 

v(t^):  specifically,  E{  v(  t^)  v 1  ( t  ^. ) }  =  R6^. 

and  P  =  state  covariance  (uncertainty)  matrix  defined 

as  the  matrix  of  products  from  P.  .  =  r.  .a. a., 

ij  ij  i  J 

r^j  defined  as  the  correlation  coefficient 
which  is  equal  to  one  for  i=j ,  and  o  defined 
as  the  standard  deviation  of  the  appropriate 
state 


With  these  defined,  the  linear  measurement  model  for  the  pro¬ 
posed  correlator/filter  tracker  can  be  defined. 

For  this  effort,  the  measurement  update  model  is  defined 
as 


l(t . )  =  HxJt.  )  +  v(t  ) 

1  —  r  l  —  1 

where 

„  [1  0  0  0  0  0  1  0' 

-  01000001 


(4-4) 


(4-5) 


Definition  of  the  model  in  this  manner  yields  the  Kalman 
filter  measurement  update  equations  in  linear  form.  Specif¬ 
ically,  these  are 


K(ti)  =  P(t7)HT[HP(tpHT  +  RU.)]"1  (4-6) 

xF(t!)  =  xF(t7)  +  K(ti)[z(ti)  -  HxF(t7)]  (4-7) 

P(t+)  =  P(tT)  -  K(t.)HP(tT)  (4-8) 

X  X  X 


where  H,  as  in  Equation  (4-4).  is  defined  to  be  time- invariant . 

Since  Equations  (4-6)  through  (4-8)  are  to  be  wholly 
implemented  within  the  optical  Kalman  filter  (0KF)  described 
in  Chapter  V,  no  further  differentiation  of  filter  measure¬ 
ment  update  operation  is  required.  However,  a  statistical 
characterization  of  v(t^)  must  be  accomplished  in  order  to 
specify  the  model  upon  which  the  update  equations  are  based 
completely . 
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Previous  work  in  this  area  (21 ,23)  has  characterised 
these  statistics  for  correlator  threshold  fractional  values 
of  0.3  and  0.5.  Since  a  threshold  fractional  value  of  0.3 
is  defined  for  this  research,  only  statistics  for  this  case 
are  reported.  Figures  26  and  27  thus  represent  previously 
obtained  (21)  error  histograms  for  the  correlation  errors 
produced  by  the  correlator  described  in  Chapter  III. 

These  error  histograms  were  obtained  in  the  following 
manner.  An  offset  between  the  target  and  template  of  0.15 
pixels  was  selected  as  being  representative  of  the  antici¬ 
pated  propagation  error  of  the  filter.  Two  sets  of  1000  runs 
each  were  executed  in  software  (one  for  the  single  hot  spot 
case  and  the  other  for  the  multiple  hot  spot  case)  to  develop 
the  correlator  statistics  for  this  offset  error.  Each  set 
used  the  truth  model  developed  in  Chapter  II  to  provide  the 
true  trajectory  information,  resulting  FLIR  intensity  pattern 
background  noises  (allowed  to  vary  over  the  full  range  of 
signal- to-noise  ratios  of  10  to  20),  and  FLIR  noises  due  to 
the  spatial  correlation  betweening  adjoining  pixels.  The 
offset  error  thus  defined  is  the  deviation  of  the  offset 
produced  by  the  correlator  from  the  correct  value  of  0.15 
pixels.  The  offset  error  produced  by  each  run  was  collected 
and  placed  in  a  bin  which  is  0.01  pixels  wide  and  within 
-0.2  pixels  to  +0.4  pixels  of  the  true  offset.  Each  histo¬ 
gram  thus  depicts  the  total  number  of  times  (out  of  1000) 
that  the  error  corresponding  to  each  bin  was  observed  for 
the  case  identified.  The  statistics  (in  pixels)  produced 
by  these  simulations  are  listed  in  Table  I. 

Table  I.  Correlation  Errors 


Use  of  the  values  contained  in  Table  I  yields  (since  the 
errors  in  each  FLIR  dimension  are  independent  of  each  other) 

*.01  759  3 

R(t.)  =  (4-9) 

0  .01  789, 

for  the  single  hot  spot  case,  and 

'.001  49  0 

R(t  )  =  ( 4- 1 C ) 

0  . 00295_ 

for  the  three  hot  spot  case.  Since,  as  outlined  earlier, 
these  statistics  were  generated  over  a  large  number  of  runs 
which  used  the  full  range  of  anticipated  conditions  for  the 
simulation,  the  resulting  matrices  of  Equations  (4-9)  and 
(4-10)  can  be  assumed  representative  for  each  type  of  target 
defined  in  Chapter  II  ?s  long  as  the  remainder  of  the  condi¬ 
tions  specified  for  the  truth  model  are  not  violated. 

Examination  of  these  statistics  reveals  three  important 
points.  First,  it  is  obvious  that  both  R(t^)  matrices  are 
required  in  the  analysis.  Thus,  the  tracker  must  be  told  the 
type  of  target  it  is  to  track  along  with  the  initial  position 
data.  Secondly,  it  is  also  easy  to  see  from  the  error  analy¬ 
sis  that  R(t^)  is  time-invariant  given  the  methodology  used 
to  develop  the  statistics.  Lastly,  the  degree  of  difference 
between  the  two  FLIR  dimension  error  uncertainties  for  the 
three  hot  spot  case  warrants  explanation.  Here,  this  effect 
is  due  to  the  fact  that  a  three  hot  spot  oar get,  as  developed 
in  Section  2.5  of  this  thesis,  would  have  an  intensity  patter 
which  is  not  radially  symmetric  about  the  center  of  the  FLIR 
image  if  it  is  orojected  so  that  the  COM  of  its  intensity 
profile  is  at  the  center  of  the  FLIR.  ‘.vith  these  three  point 
clarified,  the  entire  measurement  update  model  and  its  re¬ 
sulting  filter  equations  are  defined  for  the  targets  as  speci 
fied  by  the  truth  model  of  Chapter  II. 

The  system  dynamics  model  used  by  the  filter  must  now  be 
defined.  This  model,  outlined  in  the  final  section  of  Ohio 
chapter,  fakes  the  enhanced  state  estimate  vector  and  the 


enhanced  state  covariance  matrix,  x7(t.'  and  ?(V. 
gates  them  forward  based  on  its  knowledge  (and  th 
uncertainty  of  this  knowledge)  of  the  system  dyna 
x(t7  +  i )  and  P(t^+-j)  at  the  next  sample  period.  T 
keeps  repeating  until  weapon  operation  ceases.  T 
this  propagation  to  the  next  sample  period  has  be 
the  tracker  algorithm  itself  is  fully  defined. 

1.3  System  Dynamics  Model 

In  contrast  to  the  previously  used  first-ord 
Markov  linear  target  acceleration  model,  the  cons 
rate  target  acceleration  model  considers  the  targ 
after  projection  onto  the  FLIP,  image  plane,  to  be 
modelled  as  constant-speed,  constant  turn-rate  ma 
The  nonlinearity  in  this  model  arises  from  the  de 
of  the  target's  jerk  level  motion  (the  time  deriv 
acceleration) . 

Specifically,  the  CTP  target  acceleration  ay 

(5,3,18,29)  is 


which  implies  a  state  differential  equation 


x_(t)  =  f[x,(t),tj  +  G7Wp(  t) 


f  =  |va(t)  Vg(t)  aa(t)  ag(t)  -<^2(t)va(t)  -uj2(t) 

(-1  /t  ^)ol  (-1/ 

G^w-,  ( t )  =  jO  0  9  0  wn_(t)  wn„(t)  w,~(t)  w  (t)| 


v  x  a 


vcta3  ~  v3aq 

V..1  +  vV 


and  where 


f[xp(t),tj  =  the  nonlinear  system  dynamic 
of  the  filter  states 
cj  =  the  target's  turn-rate  on  th 
image  plane 


w^,  5  zero  mean,  white,  Gaussian  noise  composed  of 

both  filter  target  dynamics  model  uncertainties 
(subscript  DF)  and  filter  atmospheric  model 
uncertainties  (subscript  AF)  common  to  both 
FLIR  plane  dimensions  with  statistics  developed 
by  tuning 

Gp  =  dynamics  driving  noise  input  matrix 
and  2  correlation  time  assumed  by  the  filter  for  the 

atmospheric  jitter 

Note  that  this  preserves  the  original  state  matrix  definition 
expressed  in  Equation  (4-3)* 


The  statistics  of  the  noise  vector  are 


S{Wp(t)Wp(t+T) }  =£p6(T) 


(4-13) 


where  w^(t)  is  defined  as  a  four-element  column  vector  for 

— F 

G„  defined  as 

" 

.  -4x4  _ 

Definition  in  this  manner  yields  the  4x4  £p  matrix  as 


2oaf,/~af 


2oaf/taf 


where 


°DF  “  the  assumed  target  jerk  level  motion  white  pro¬ 
cess  power  spectral  density  value 

and  2  the  assumed  atmospheric  jitter  Markov-1  process 

A  r 

variance 


From  the  model,  the  equations  necessary  to  propagate  the 
state  estimates  and  state  estimate  uncertainties  forward  in 
time  (to  the  next  sample  period)  are  developed.  In  discrete 


time,  the  general  form  of  the  state  estimate  propagation  is 

Xp(tT+1)  =  xF(tt)  +  J^i+I  f  [xF(  1 1  ti)  ,  t]  dt  (4-14) 

1 

where  here  the  nonlinear  system  dynamics  function  is  condi¬ 
tioned  on  time  .  If  a  first-order  Euler  integral  approxi¬ 
mation  to  Equation  (4-14)  is  used  over  the  interval  4t=t^+^-ts 
the  expression 

xF  ( t  T  + 1 )  =  x^Ct^)  +  f  [x-,(  tt)  .tjat  (4-15) 


can  be  used  for  a  At  which  is  sufficiently  small  when  compared 
to  the  natural  transient  times  of  the  physical  system.  As 
such  is  the  case  here,  Equation  (4-15)  is  totally  computed  in 
software  to  produce  the  state  estimate  propagation.  This 
will  be  further  discussed  in  Chapter  V  where  the  method  of 
filter  operation  is  discussed  in  terms  of  efficiency. 

Similarly,  the  covariance  propagation  (state  estimate 
uncertainty  propagation)  must  also  be  written  in  terms  of  the 
most  recent  state  estimates  since 


F(ti)  =  df[x..t]/3x 


(<) 


(4-16) 


is  no  longer  known  a  priori  as  it  would  be  in  a  linear  filter. 
The  definition  of  F(t^)  in  Equation  (4-16)  as  a  piecewise  con¬ 
stant  function  again  is  valid  due  to  the  relative  rates  of 
variation  of  F(t.)  with  At. 

The  standard  linear  covariance  propagation  equation  can 
thus  be  used  with  the  restriction  that  the  parameters  of  in¬ 
terest  be  recalculated  each  sample  period.  Specifically, 
this  equation  is 


—  (  ^  i+ 1  ^  ~  — F^'i+1  fti^  — ^i^— F^i+1  ’  ^i^  +  ^FD^i+1  ’  *i^  (4-1  <) 
where  the  variables  of  interest  are  developed  as  follows. 
Using  the  same  first-order  Euler  approximation  as  in  the 
state  estimation  propagation  development, 


$F(tn. +  1,ti)  =  expfF(ti)At’J  (4- IS) 

For  the  nonlinear  system  dynamics  function  specifically  de¬ 
fined  by  Equation  (4-12),  it  can  readily  be  seen  tha  +  calcu¬ 
lation  of  ( t..  +  ^  ,  t.  )  using  the  strict  interpretation  of 
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Equation  (4-18)  would  impose  a  large  computational  burden  on 
the  system.  This  is  particularly  true  since  it  is  to  be  done 
in  software.  However,  a  careful  analysis  of  the  matrices  of 
Equation  (4-17)  reveals  that  any  terms  not  on  the  diagonal  of 
$F  will  have  zero  contribution  to  P.  This  is  because  p  itself 
is  diagonal  due  to  the  independence  of  all  of  the  system  un¬ 
certainties  with  each  other.  With  this  knowledge,  the  upper 
left  6x6  portion  of  the  matrix  is  truncated  to  first-order 
terms,  resulting  in  the  expression 

£F(t.+1,t.)  =  I  +  F(t.)At  (4-19) 

The  balance  of  F(t^),  being  associated  with  the  two  atmo¬ 
spheric  states,  is  time-invariant  and  can  be  determined  in 
its  exact  closed  form.  Appendix  B  provides  a  full  analysis 
of  this  development. 

For  the  same  reasons,  the  £FF  matrix  is  also  simplified. 
Formally,  its  definition  is 

Wh+I’h*  =  Jttiit1  iF(tit1"T)£F2F<T>£F£F(ti+1'T)  dT  (4~20) 

which  must  be  recalculated  every  sample  period  since  the  F 
upon  which  6  is  based  is  not  known  a  priori.  However,  since 
the  time  interval  At  is  short  compared  to  the  system  tran¬ 
sients,  previous  research  has  shown  (8,29)  that  the  defini- 

tl0n  %D  =  Gp^At  (4-21) 

is  valid  for  the  filter  specified.  Given  the  structures 
already  identified  for  both  Q_  and  G„,  it  is  obvious  that 
once  the  values  are  specified  (through  tuning)  for  £F , 
becomes  constant  with  values  of  C^At, 

x 

Processing  of  the  filter’s  propagation  cycle  is  thus 
accomplished  by  using  software  to  execute  Equations  (4-16) 
and  (4-19)  on  an  iterative  basis.  Software  will  also  be  used 
to  calculate  all  constant  values  initially  upon  each  operation 
of  the  tracker.  The  OKF  will  then  implement  Equation  (4-17) 
to  provide  the  state  estimate  covariance  propagation. 
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With  these  developments,  the  entire  tracker  has  been 
mathematically  specified.  The  remaining  task  of  this  thesis 
prior  to  analyzing  the  tracker's  performance  is  the  complete 
development  of  the  optical  Kalman  filter.  The  next  chapter 
is  devoted  to  this  requirement,  as  well  as  tieing  the  entire 
system  together  in  a  cohesive  form  prior  to  its  analysis. 


V.  Optical  Kalman  Filter 


5.1  Introduction 

With  the  filter  models  fully  developed  and  analysed,  the 
method  of  implementation  of  the  filter  algorithms  must  be 
completely  described  and  characterized.  This  requirement  is 
particularly  true  given  the  optical  nature  of  the  proposed 
solution.  The  processor  analysis  must  consider  the  accuracy 
of  the  results  generated  by  the  optical  bed  versus  those 
produced  by  a  software  implementation  (given  the  accuracy  re¬ 
quired  by  the  problem  dynamic  range)  in  addition  to  system 
throughput  considerations ,  Successful  demonstration  of  spe¬ 
cification  compliance  will  thus  constitute  satisfaction  of 
the  major  basis  of  this  thesis. 

Several  optical  architectures  were  considered  as  candi¬ 
dates  for  the  optical  processor.  However,  the  tradeoff  anal¬ 
yses  outlined  in  Chapter  I  which  resulted  in  the  choice  of 
the  pipelined  iterative  optical  systolic  array  architecture 
are  negligible  compared  to  both  a  description  of  the  archi¬ 
tecture  itself  and  the  way  the  architecture  implements  the 
operations  required  to  perform  the  filter  algorithms.  For 
this  reason,  the  discussions  of  this  chapter  are  limited  to 
these  areas  and  the  required  analyses  thus  generated. 

To  accomplish  these  developments  in  the  most  logical 
manner,  the  initial  thrust  of  this  chapter  (Section  5.2)  will 
be  limited  to  development  of  the  optical  architecture  and  its 
method  of  implementing  the  operations  required  by  the  general 
Kalman  filter  algorithm.  To  maintain  a  level  of  presentation 
sufficiently  basic  to  allow  understanding  by  engineers  unfa¬ 
miliar  with  optics,  device  physics  will  not  be  fully  developed 
within  the  section  but  will  be  referenced  for  optional  study. 
Experienced  optical  engineers  may  wish  to  substitute  Reference 
2,  upon  which  this  architecture  is  based,  in  lieu  of  this 
section . 


The  final  section  of  this  chapter  contains  all  th 
yses  required  to  demonstrate  processor  specification  c 
ance  after  all  of  the  specific  real-world  components  h 
been  identified  and  characterized.  Sufficient  analyse 
provided  to  present  a  clear  comparison  between  the  pro 
optical  approach  and  a  software  implementation  of  the 
algorithms  developed  in  Chapter  IV. 

ipelined  Iterative  Optical  Svstolic  Arrav  Seveio 


The  decision  to  employ  an  optica^,  processor  as  the  solu¬ 
tion  of  choice  is  chiefly  due  to  the  inherent  ability  of  op¬ 
tical  components  to  handle  a  large  degree  of  data  in  parallel 
without  requiring  component  duplication.  The  particular 
architecture  selected,  the  pipelined  iterative  optical  sys¬ 
tolic  array,  fully  exploits  this  capability  as  well  as  ex¬ 
ploiting  the  advantages  presented  by  data  pipelining  and  the 
iterative  nature  of  the  process  itself. 

Prior  to  developing  the  architecture,  clarity  requires 
that  each  optical  device  type  used  within  the  architecture 
be  fully  characterized  in  order  to  illustrate  clearly  the 
data  flow  through  the  assembled  architecture.  Since  the 
architecture  itself  is  extremely  simple,  only  four  optical 
devices  require  characterization:  light  sources,  acousto-op¬ 
tic  light  modulators  (AOMs),  simple  lenses,  and  optical  detec¬ 
tors.  The  description  of  each  particular  device's  character¬ 
istics  will  be  presented  in  the  order  of  their  occurrence  in 
the  architecture  and  will  be  tied  back  to  the  previous  device 
characteristic  descriptions  so  as  to  maintain  a  logical  data 
flow  which  is  easily  relatable  to  the  overall  architecture. 

The  first  optical  component  type  required  by  the  archi¬ 
tecture  is  a  linear  array  of  light  sources  in  order  to  pro¬ 
duce  a  parallel  signal  representation  in  the  optical  domain. 
This  linear  array  requirement,  as  well  as  the  size  constraints 
imposed  by  the  modulator  (as  will  be  seen),  limits  the  choice 
of  optical  components  to  two  distinct  but  related  semicon¬ 
ductor  device  types  (25,28,30),  either  light  emitting  diodes 


(LEDs)  or  laser  diodes  (LDs)  ,  combined  as  discrete  units  or 
a  prepackaged  array.  While  this  limitation  may  appear  to  be 
severe,  the  range  of  options  available  to  the  system  archi¬ 
tect  is  actually  quite  large. 

LEDs  are  classified  as  incoherent  sources  of  light  and 
are  available  in  a  wide  range  of  maximum  average  power  levels, 
continuous  wave  and  pulsed  operational  modes,  and  output 
bandwidths  (light  spectral  content,  pulse  repetition  rate, 
and  pulse  width).  However,  their  incoherency,  an  output 
distributed  over  some  range  of  light  frequencies,  produces 
five  inherent  device  limitations.  These  limitations  are 
directly  related  to  the  light  frequency  range  of  the  device 
and  can  limit  their  applicability  in  optical  processing 
architectures.  These  limitations  include  constraints  on 
(1)  the  maximum  pulse  repetition  rate  and  (2)  the  minimum 
pulse  width  of  the  device.  The  component  density  is  also 
limited  since  (3)  an  LED's  projected  field  of  view  is  much 
greater  than  that  of  a  coherent  source.  Additionally,  these 
devices  are  highly  susceptible  to  (4)  spatial  and  (5)  spectral 
output  variations  over  time  due  to  thermal  effects,  although 
these  drawbacks  can  be  somewhat  mitigated  by  cooling.  Despite 
the  fact  that  these  limitations  do  sometimes  prevent  the  use 
of  LEDs  in  optical  architectures,  use  of  these  devices  is  as 
often  mandated  (although  not  in  this  case)  by  their  advantage 
in  available  maximum  average  output  power  over  other  semicon¬ 
ductor  source  types. 

Laser  diodes,  on  the  other  hand,  are  coherent  devices 
since  they  employ  an  actual  laser  cavity  to  produce  light. 
Cavity  losses,  principally  the  extinction  of  frequencies 
not  harmonic  to  the  cavity,  thus  limit  their  maximum  average 
output  power,  particularly  when  compared  to  LEDs  having 
identical  active  regions  in  both  size  and  material.  Although 
this  limitation  sometimes  is  a  serious  drawback  to  their  use 
in  optical  architectures,  the  fact  that  the  coherency  of  a 
LD  allows  extremely  narrow  pulse  widths  resulting  in  high 
peak  power  outputs  permits  their  use  in  the  majority  of 
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non-integrating  architectures.  Additional  benefits  as- 
cribable  to  coherency  include  greater  pulse  repetition  rates, 
greater  device  density  due  to  the  much  narrower  pro  jeered 
FOV,  and  minimization  of  thermal  effects  due  to  their  ability 
to  maintain  a  coherent  output.  For  these  reasons,  LDs  are 
often  the  only  sources  available  for  high-speed  optical  pro¬ 
cessing  applications. 

Either  of  these  device  types,  when  configured  in  a  lin¬ 
ear  array,  formulate  a  parallel  representation  in  a  spatial 
sense  of  either  a  group  of  analog  values  or  a  digital  bit 
stream.  Since  the  projected  paths  of  the  signals  formulated 
by  these  devices  are  parallel  or  at  least  separated  over  the 
region  of  interest,  as  illustrated  in  Figure  28,  modulation 
of  the  values  represented  by  each  light  source  can  also  occur 
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Figure  28  - 


Parallel  Spatial  Sense  Signal  Representation 
in  (a)  LEDs  (b)  LDs 


in  parallel.  For  the  architecture  chosen,  this  modulation, 
as  noted  in  Chapter  I,  is  accomplished  through  use  of  a 
single  AOM  ( 30) . 

The  first  requirement  for  modulating  a  spatially  parallel 
signal  representation  is  that  the  signal  loaded  into  the  mod¬ 
ulator  be  spatially  matched  to  the  input  signal  produced  by 
the  LED/LD  array.  In  an  AOM,  this  is  accomplished  by  serially 
injecting  a  series  of  bulk  acoustic  waves  (discrete  frequency 
packets)  into  a  crystal  which  has  been  carefully  oriented 
along  a  specific  crystal  axis.  These  injected  discrete  fre¬ 
quency  packets  are  often  separated  by  nulls,  as  illustrated 
in  Figure  29,  to  preclude  crosstalk  between  values,  and  have 
power  levels  corresponding  to  each  particular  value.  But 

LED/LD  AOM  AOM  Cell  Position 


Figure  29  -  Matched  Parallel  Signal  Representation 

AOMs  inherently  possess  two  serious  limitations  which  must  be 
overcome.  These  are  acoustic  attenuation  within  the  crystal 
and  the  acoustic  beam  spread  as  a  function  of  propagation. 

These  two  limitations  present  a  number  of  tradeoff  con¬ 
siderations  to  the  architect.  However,  the  limitation  im¬ 
posed  by  acoustic  attenuation  can  be  overcome  by  using  a 


source  (LD  or  LED)  which  has  sufficient  output  power  both  to 
represent  the  signal  input  to  it  accurately  and  to  compensate 
for  the  attenuation  of  the  signal  in  the  ACM.  This,  of 
course,  presupposes  that  the  system  throughput  rate  combined 
with  the  velocity  of  sound  (v  )  in  the  ACM  produces  a  work¬ 
ing  aperture  (spatial  extent  of  the  parallel  representation) 
which  is  of  sufficient  size  to  accommodate  the  required  source 
devices.  Additionally,  the  choice  of  crystal  and  crystal  axis 
determines  to  some  extent  the  degree  of  correction  obtainable 
to  counter  the  acoustic  beam  spread  limitation.  This  is  due 
to  the  fact  that  choice  of  transducer  geometry  only  collimates 
(forms  a  spatially  parallel  representation  of)  the  propagation 
vectors  of  the  acoustic  waves  at  the  AOM  input  and  does  not 
affect  distortion  caused  by  the  sources  identified  above. 

Assuming  for  the  moment  (until  satisfactory  demonstra¬ 
tion  in  Section  5.3)  that  these  limitations  have  been  over¬ 
come  so  that  the  AOM  output  response  is  flat  across  the  en¬ 
tire  working  aperture  (i.e.,  the  response  to  a  unit  input  to 
both  the  AOM  and  each  element  of  the  LED/LD  array  is  uniform 
across  the  entire  spatial  representation  of  the  AOM  output), 
the  process  of  analog  multiplication  is  described.  Recall 
that  the  bulk  acoustic  waves  are  amplitude  modulated  (in  terms 
of  acoustic  power  per  unit  area)  discrete  frequency  packets 
which  propagate  along  the  longitudinal  axis  of  the  AOM  as 
illustrated  in  Figure  29.  Then  for  any  specified  acoustic 
bandwidth  each  discrete  frequency  packet  can  be  composed  of 
a  superposition  of  several  amplitude  modulated  frequencies. 

Each  acoustic  frequency  (w„)  present  within  the  AOM  can  then 
interact  with  the  input  light  frequency  (uj)  in  the  following 
manner . 

Restricting  for  the  moment  the  modulating  frequency  to 
the  AOM  acoustic  center  frequency,  a  Bragg  diffraction  condi¬ 
tion  (30)  is  purposely  implemented  as  illustrated  in  Figure 
30.  In  terms  of  the  wavelengths  involved,  this  Bragg  condition 
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Figure  30  -  AOM  Bragg  Diffraction  Condition 

exists  when  the  path  length  AO+OB  is  equivalent  to  the  optical 
wavelength  within  the  crystal  (the  free  space  wavelength.  A, 
divided  by  the  crystal's  index  of  refraction,  n)  .  Since 


A  sin  3  =  A/2n 
s 


(5-1) 


where  A^  is  the  wavelength  of  the  acoustic  wave  within  the 
s 

crystal,  the  precomputable  Bragg  angle  required  to  set  up 
the  diffraction  is 


3  =  sin 


.  -1  A 

m  rt — 7“  =  31. 

!_2nA  J 


1  I  Av, 


since  Ag  is  defined  as  the  velocity  of  sound  in  the  crystal 
(vo)  per  unit  sound  frequency  (v  ).  In  terms  of  frequency, 
the  diffracted  beam's  angular  propagation  (vector)  character¬ 
istic  is  now  described  by  to+ui  .  Relaxing  the  fixed  frequency 
assumption,  the  relationship  [with  the  input  Bragg  angle  as 
defined  in  Equation  (5-2)  maintained]  reduces  to 

A0  =  — -Av  (5-3) 

nv  s 
s 

where  Av  is  the  difference  between  the  input  sound  frequency 
and  the  acoustic  center  frequency  and  A0  is  the  resulting 
angular  output  difference  from  the  Bragg  angle  (30).  Noting 
that  Av^  can  be  either  positive  or  negative,  each  discrete 
sound  frequency  packet  input  to  the  AOM  spatially  modulates  t: 


1 


incident  light  over  an  angular  range  defined  by 
acoustic  bandwidth.  This  is  expressabie  as 
I, 


tee 
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diffracted  light  exitance  (emitted  power 
unit  area) 

incident  light  irradiance  (power  per  unit 
area  onto  a  surface) 

diffraction  efficiency  of  the  AOM  or  the 
Bragg  efficiency 
wavelength  of  light  in  uni 
the  acoustic  interaction  length  (generall 
taken  to  be  the  width  of  the  acoustic  bea 
through  which  light  would  pass  if  it  were 
input  perpendicular  to  the  direction  of 
acoustic  propagation) 

the  diffraction  figure  of  merit  of  the 
crystal  relative  to  water 


and  -^A^oustic  ~  acoustic  intensity  within  the  crystal 


(acoustic  power  per  mm  area) 


As  the  input  light  intensities  are  also  amplitude  modulate 
it  is  easy  to  see  that  for  the  analog  case  it  is  the  indiv 
dual  products  of  the  values  represented  by  the  incident  li 
power  and  each  of  the  frequencies  present  in  the  discrete 
frequency  packet  that  are  angularly  separated,  whereas  for 
the  digital  case,  multiple  two-bit  AND  functions  are  pro¬ 
duced.  cince  this  occurs  at  every  modulation  point  in  the 
AOM  parallel  representation ,  discrimination  between  proauc 
is  impossible  without  further  spatial  modulation.  In  the 
chosen  architecture,  this  required  spatial  modulation  is 
accomplished  by  a  . •  i m tie  lens. 
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Thus,  for  the  case  of  a  lens  oriented  so  that  its  diameter  is 
perpendicular  to  the  propagation  vector  (from  an  AOM)  re¬ 
sulting  from  the  Bragg  angle  and  having  all  products  modulated 
using  either  positive  or  negative  multiples  of  a  fixed  fre¬ 
quency  Av  about  the  AOM  acoustic  center  frequency,  the  sep- 
aration  between  any  two  product  sums  is  expressible  as 

Ax  =  f  tan  A0  (5-7a) 

- f  tan[^]avs  (5-7b) 

where  the  multiples  of  Av  over  the  product  sum  range  are 

s 

only  constrained  by  the  AOM  acoustic  bandwidth. 

There  are,  however,  some  constraints  for  the  utilization 
of  these  relationships.  Since  they  are  based  on  thin  lens 
formulae,  the  lens  diameter  must  be  oriented  perpendicular 
to  the  Bragg  angle  0  to  produce  the  A0  relationship  expressed 
in  Equation  (5-7)  accurately.  Additionally,  due  to  lens 
aberration  properties,  minimization  of  error  is  achieved  only 
when  the  inner  portion  (roughly  corresponding  to  the  semi¬ 
circle  described  by  half  the  lens  diameter)  of  the  lens  is 
used  and  the  lens  F-number  (lens  focal  length  divided  by  the 
lens  diameter)  is  relatively  high  (>_  f/4)  to  preclude  the 
effects  of  large  angles.  Moreover,  all  of  these  constraints 
become  increasingly  more  difficult  to  satisfy  as  the  system 
throughput  rate  and/or  array  size  is  increased. 

Now  that  the  specific  product  sums  have  been  focused  at 
the  same  point  in  space,  a  linear  detector  array  (25 ,28,30) 
is  required  to  convert  the  optical  signal  back  into  an 
electronic  signal  accurately.  Due  to  size  constraints,  the 
choice  is  limited  to  semiconductor  devices,  where  the 
parameters  of  interest  are  the  device  quantum  efficiency  for 
the  particular  light  wavelength  used,  the  device  minimum 
detectable  power,  and  the  maximum  (or  saturation)  power  of 
the  device.  Since  the  final  two  parameters  are  directly 
related  to  the  size  of  the  device's  active  area,  the  major 
constraint  is  once  again  maintaining  sufficient  lateral  sep¬ 
aration  to  accomodate  component  size. 
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With  the  devices  incorporated  within  the  architecture 
defined,  a  step-by-step  development  of  the  linear  algebra 
operations  required  for  a  Kalman  filter  algorithm  must  be 
accomplished.  Although  only  illustrated  here  for  square 
matrices  and  vectors  of  dimension  three,  the  architecture 
and  processing  steps  are  directly  expandable  to  any  higher 
dimension  within  the  aforementioned  component  constraints. 

The  least  complex  linear  algebra  required  for  a  Kalman 
filter  is  the  formation  of  a  matrix-vector  product,  ex¬ 
pressible  as 


aii  ai2  ai3 

'  bi  ' 

’  Cl 

a2  1  a2  2  a2  3 

b2 

= 

c2 

a 3  i  a 3 2  a3  3 

.  b3  - 

-  C3  _ 

Using  the  architecture  as  illustrated  in  Figure  32  (with 
the  Bragg  angle  not  depicted  for  clarity)  and  superimposing 
the  matrix  column  values  in  the  sequence  outlined  in  Table 
II,  the  complete  matrix-vector  product  is  available  in  (2N-1) 
time  intervals  for  a  matrix  of  dimension  MxN  and  a  vector  of 
dimension  N.  For  tabular  purposes,  AOM  modulation  frequen¬ 
cies  are  expressed  as  diffraction  orders  where  each  full  dif¬ 
fraction  order  is  defined  to  be  a  multiple  of  Av  about  the 

s 

AOM  acoustic  center  frequency.  Thus  for  an  odd  number  of 
superimposed  frequencies,  the  diffraction  orders  are  expressed 
as 

0,  ±1,  ±2,  ±3,  ...,  ±  - 

where  S  is  the  number  of  superimposed  frequencies,  whereas 
for  the  case  of  an  even  number  of  superimposed  frequencies, 
the  diffraction  orders  are 

0,  ±0.5,  ±1.5,  ±2.5,  . . .  ,± 


Use  of  this  technique  maintains  a  symmetric  diffraction 
pattern  about  the  AOM  acoustic  center  frequency  and  thus 
minimizes  potential  error  sources  as  well  as  maintaining  the 


Table  II.  Matrix-Vector  Formation  Sequence 


Time  Actions 


Outputs 


an  at  +1,  a2 1  at  0,  and  a33  at  -1 
input  to  AOM  at  Ti 


T2  1st  column  of  A  propagates  to  x2 
2nd  column  of  A  input  to  AOM  §  ti 
using  same  superposition  frequencies 


1 st  column  of  A  propagates  to  t3 
2nd  column  of  A  propagates  to  t2 
3rd  column  of  A  input  to  AOM  @  Ti 
using  same  superpostion  frequencies 


1 st  column  of  A  propagates  to  x4 
2nd  column  of  A  propagates  to  x3 
3rd  column  of  A  propagates  to  x2 


1 3 t  column  of  A  propagates  to  x5 

2nd  column  of  A  propagates  to  x<* 

3rd  column  of  A  propagates  to  x3 

bi  input  to  LED/LD  1 

b2  input  to  LED/LD  2 

b3  input  to  LED/LD  3 


-  •>  •  •  V*  •  •  -  * 


,cvj 


relationship  expressed  in  Equation  ( 5 — 7b > .  Therefore,  in 

Table  II  (and  other  appropriate  tables  in  this  section) , 

"an  at  +1"  means  "an  at  diffraction  order  +1"  or  an  at 

a  frequency  of  v  +  (1)Av  ",  and  so  forth. 
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Note:  Drawing  Not  to  Scale 
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Figure  32  -  Matrix-Vector  Multiplication  (at  t=T5) 


Detectors 


Although  it  can  be  readily  seen  that  a  matrix-vector 
product  could  be  available  at  either  T3,  T4»  or  T5»  use  of 
the  same  architecture  for  matrix-matrix  or  higher  order 
operations  requires  the  timing  illustrated.  This  can  be 
demonstrated  by  merely  tabulating  the  required  operations 
for  matrix-matrix  multiplication  (here  A B  = C)  by  the  same 
architecture.  This  process,  illustrated  by  the  combination 
of  Table  III  and  Figure  33a-c,  again  requires  only  (2N-1) 
time  intervals  to  produce  the  desired  outputs.  In  general, 
the  multiplication  of  two  matrices  of  dimensions  (MxN)  and 
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Table  III  -  Matrix-Matrix  Multiplication  Sequence 


Time  Actions 


1 st  column  of  A  input  to  AOM  @  Ti 
using  same  superposition  frequencies 
as  matrix-vector  case 

1 st  column  of  A  propagates  to  t z 
2nd  column  of  A  input  to  AOM  @  Ti 
using  same  superposition  frequencies 

1 st  column  of  A  propagates  to  ta 
2nd  column  of  A  propagates  to 
3rd  column  of  A  input  to  AOM  @ 
using  same  superposition  frequencies 
bn  input  to  LED/LD  3 
bz i  input  to  LED/LD  4 
b3i  input  to  LED/LD  5 

as  as  -  j  SB  —  a  as  ^b  — —  ss  a  ^a  s  a  si 

1 3t  column  of  A  propagates  to  T4 

2nd  column  of  A  propagates  to  T3 

3rd  column  of  A  propagates  to  xz 

bi2  input  to  LED/LD  2 
b22  input  to  LED/LD  3 
b 3  2  input  to  LED/LD  4 


1 st  column  of  A  propagates  to  ts 

2nd  column  of  A  propagates  to  x4 

3rd  column  of  A  propagates  to  13 

bi3  input  to  LED/LD  1 

b2  3  input  to  LED/LD  2 

bss  input  to  LED/LD  3 


Outputs 


Note:  Drawings  Not  to  Scale 


FT 


Figure  33  -  Snapshots  of  Matrix-Matrix  Product  Formation 

at  (a)  T3  (b)  T4  (c)  T5 


(NxM)  requires  M  frequencies,  M  detectors,  (2N-1)  LEDs/LDs, 
and  (2N-1)  time  intervals  of  which.  (N-1)  time  intervals  is 
dead  or  delay  time. 

If  the  linear  algebra  operations  required  by  the  Kalman 
filter  algorithm  were  limited  solely  to  the  product  formations 
already  illustrated,  architectural  development  could  stop  at 
this  point.  However,  the  necessity  of  being  able  to  perform 
matrix-matrix-matrix  multiplies  and  matrix  inversions  in  an 
efficient  manner  requires  the  addition  of  three  externally 
engageable  electronic  devices  to  the  optics  already  illus¬ 
trated.  The  efficiency  requirement  precludes  repetition 
of  matrix-matrix  operations  to  accomplish  these  higher  order 
operations  if  better  alternatives  are  available. 

For  example,  the  architecture  required  to  support 
matrix-matrix-matrix  product  formation  is  illustrated  in 
Figure  34  and  requires  inclusion  of  a  sample  and  hold  (S/H) 
device  in  the  added  feedback  loop.  Noting  that  the  optical 
bed  is  identical  to  those  discussed  previously,  only  a 
tabular  summation  of  processing  steps  for  the  product 
E  =  ABD  (where  kQ  =  C  and  E  =  £D)  is  presented  in  Table  IV. 
From  Table  IV,  it  is  obvious  that  the  overall  computation 
time  required  for  the  complete  formation  of  E  is  (3N-1) 
time  intervals  of  which  (N-1)  time  intervals  is  again  dead 
or  delay  time.  Further  matrix  multiplications  are  also 
achievable  using  the  same  algorithm  with  each  additional 
matrix  requiring  an  additional  N  time  intervals  to  form  the 
complete  product. 

Rather  than  illustrating  further  pipelining  beyond 
matrix-matrix-matrix  multiplication,  matrix  inversion  using 
the  same  optical  package  is  outlined.  Here,  a  modified 
Richardson  algorithm  (2)  to  solve  C  =  H  B  for  B  = 
without  explicitly  computing  the  matrix  inversion  H~  1  is 
utilized.  Use  of  this  algorithm  exploits  the  pipelined  itera¬ 
tive  nature  of  the  algorithm  as  follows.  Within  this  concept, 


Figure  34  -  Architecture  for  Matrix-Matrix-Matrix  Multiplier 


Table  IV  -  Time  History  of  Figure  34  Components 
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Detector  Outputs 
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the  iterative  algorithm 


■2k+-i  -  ji  -  wii]3k  +  2 


(5-9) 


is  used  where  u>  is  the  acceleration  parameter  that  is  chosen 
to  speed  convergence  of  the  algorithm.  When  the  results 

on  two  successive  iterations  k  and  k+1  are  equal. 
Equation  (5-9)  reduces  to  the  solution 


B  =  H~ 1 C 


(5-10) 


The  full  architecture,  illustrated  in  Figure  35,  realizes 
the  algorithm  by  modifying  Equation  (5-9)  into 

Bk+1/u>  =  [i/co  -  HjBk  +  C  (5-11) 

where  [i/to  -  H]  is  electronically  precalculated  for  efficiency 
Realizing  that  H  is  known  and  fixed  and  that  I/oj  is  easily 
computable  by  a  simple  scaling  of  its  elements,  £l/w  -  Hj 
is  thus  written  as  a  matrix  A  which  is  assumed  known  and 
fixed  for  a  given  value  of  to,  yielding 


Skt1  ■  C]  (5-12) 

The  optical  system  thus  performs  the  matrix-matrix  multipli¬ 
cation  A  Bk,  the  matrix  G  is  added  one  row  at  a  time  in  real 
time  (in  a  resistive  adder)  to  form  [A  Bk  +  C~j  ,  and  this  matrix 
summation  is  then  electronically  multiplied  by  to  as  shown  in 
the  figure  to  form  the  next  Bk+^  iterative  input  to  the  AOM . 
Table  V  clearly  shows  the  time  history  of  the  algorithm,  where 
for  notational  simplicity  the  initial  k=0  matrix  B0  is  denoted 
by  elements  b„  .  Bi  by  b '  ,  and  B2  by  b' ’  .  As  can  be  seen 
in  the  table,  the  initial  b  inputs  to  the  AOM  originate 
from  an  initial  estimate  described  in  Reference  2  and  there¬ 
after  all  future  AOM  inputs  come  from  prior  iterations.  This 
again  yields  an  ideal  pipelining  of  operations  as  there  is 
no  dead  or  delay  time  after  the  initial  (N-1)  AOM  cell  delay. 
Given  the  filter  matrices  defined  in  Chapter  IV,  the  only 
inversion  required  for  the  problem  of  this  thesis  [Equation 
(4-6)]  is  only  of  dimension  (2x2).  By  defining  C  =  I_2 , 
the  B0  matrix  and  acceleration  parameter  to  are  easily 
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Figure  35  -  General  OKF  Architecture  (Matrix  Inversion  Case 
°  Detailed) 

.  Table  V  -  Time  History  for  Matrix  Inversion 


T1  T2  T3  T4  T5  T6  T7  T8  Parameter 


Comments 


Detector  Outputs 
D  =  AB 
D'  =  AB' 


Matrix  £ 
Inputs  to 
the  Adder 


B'=  3k+1  and 
B  ,T~  Bk+2  Inputs 

to  the  S/H 


Freq-Muxed 
Inputs  =  AOM 
Incuts  from  T4  on 


m 


defined  (2)  to  produce  a  usable  inversion  product  (where  the 
two  iteration  values  are  equivalent  within  the  processor 
accuracy)  after  only  (4N-1)  cell  operations.  This  specifi¬ 
cation  thus  completes  the  time  history  definitions  for 
all  of  the  required  Kalman  filter  operations. 

It  is  evident  that  the  general  architecture  illustrated 
in  Figure  35  can  also  accomplish  all  the  linear  algebra 
required  for  Kalman  filter  operations.  It  has  thus  been 
chosen  as  the  basic  architecture  for  processing  the  filter 
algorithms  identified  in  Chapter  IV.  However,  prior  to 
individual  component  selection  and  identif ication ,  the  method 
by  which  the  processor  handles  bipolar-valued  data  requires 
discussion.  This  requirement  arises  from  the  inability  to 
detect  the  potential  fields  due  to  light.  Since  only  power 
(or  power  per  unit  area)  can  be  detected  in  the  optical 
domain,  there  is  thus  no  capability  to  represent  negative 
numbers  directly  without  a  bipolar  algorithm. 

Since  the  overriding  concerns  in  any  bipolar-valued  data 
handling  algorithm  are  the  pipelining  of  data  and  operations 
and  the  system  throughput  rate,  utilization  of  a  feedback 
loop  requiring  extensive  A/D  and  D/A  conversions  must  be 
avoided  if  at  all  possible.  The  simplest  method  of  accom¬ 
plishing  this  algorithm  within  the  constraint  thus  becomes 
performing  the  necessary  matrix  algebra  in  the  analog  system 
to  sufficient  accuracy  as  the  rest  of  the  processor. 

Consider  first  the  multiplication  of  two  bipolar  scalars 
a  and  b.  Allowing  each  scalar  to  be  represented  in  bipolar 
i urm  using  a  ,  a  ,  b+ ,  and  b"  notations,  a  positive  scalar  is 
expressed  using  the  notation  whereas  a  negative  scalar  is 
denoted  by  the  ~  notation.  This  form  thus  permits  a  repre¬ 
sentation  of  both  positive  and  negative  scalars  which  is 
always  positive  by  defining  the  unused  half  of  the  bipolar 
representation  to  be  zero.  This  yields  the  required  bipolar 
multiplication 

ab  =  (a+  -  a  ) ( b+  -  b  )  =  (a+b+  +  a~b~)  -  (a~b’  +  ab”)  (5-13) 
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where  the  positive  and  negative  parts  of  the  product  ab  are 

(ab)+  =  (a+b++  a-b-)  (5-l4a) 

(ab)-  =  (a-b++  a+b-)  ( 5  -  1 4  b ) 

and  the  bipolar  output  product  is 

ab  =  ( ab) +  -  (ab)-  (5-15) 

It  is  evident  that  one  of  the  two  terms  of  Equation  (5-15) 
will  always  be  zero.  Extending  this  technique  to  the  case 
of  matrix-matrix  multiplication  (here,  A  B  =  _C)  where  the 
matrix  elements  are  bipolar  valued,  each  element  a  of  the 
matrix  is  arranged  as  a  (2x2)  submatrix  and  each  element  of 
is  arranged  as  a  two-element  column  vector.  For  example,  in 
the  case  of  N=M=2, 
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where  each  element  c  of  the  output  matrix  is  likewise  repre- 

mn  r 

sented  as  a  two-element  column  vector  as  shown.  As  before, 
one  of  the  two  elements  of  each  column  vector  of  _C  will  always 
be  zero  and  all  input  and  output  elements  will  always  be 
positive  or  zero. 

This  algorithm  lends  itself  to  direct  pipelining  and 
incorporation  into  all  the  linear  algebra  processing  algorithms 
since  the  matrix  output  _C  is  in  the  required  form  of  B  and 
can  thus  be  fed  back  to  the  processor.  The  algorithm  execu¬ 
tion  sequence  is  altered  as  illustrated  in  Table  VI,  resulting 
in  an  additional  N  operations  of  both  processing  and  dead  or 
delay  time  over  each  non-bipolar-valued  case.  The  net  effect 
of  this  algorithm  on  the  architecture  is  thus  merely  an 
increase  in  the  processor  size:  (3N-1)  LEDs/LDs,  2M  fre¬ 
quencies,  and  2M  detectors  are  now  required.  However,  its 
direct  integration  into  the  flow  and  pipelining  of  data  and 
operations  yields  significant  increases  in  the  potential 
sytem  throughput  rate. 
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Table  VI.  Bipolar-Valued  Matrix-Matrix  Sequence 


Time 

Actions  j 

T1 

+ 

an 

§  +1.5.  a"u  §  + .  5 ,  a2  i 

§  -.5, 

and  a"2 1  §  -1.5  input  to  AOM  at  t1 

T2 

1  st 

column  of  A  propagates 

to  t2 

2nd  column  of  A  input  to  AOM  at  t1 
same  superposition  frequencies 

using 

T3 

1  st 

column  of  A  propagates 

II 

II 

H 

II 

O 

II  -P 

2nd 

column  of  A  propagates 

to  x2 

3rd  column  of  A  input  to  AOM  at  t1 
same  superposition  frequencies 

using 

T4 

1st 

column  of  A  propagates 

to  x4 

2nd 

column  of  A  propagates 

to  t3 

3rd 

column  of  A  propagates 

to  t2 

4th  column  of  A  input  to  AOM  at  t1 
same  superposition  frequencies 

using 

bn 

input  to  LED/LD  2 

bn 

input  to  LED/LD  3 

b+2  l 

input  to  LED/LD  4 

b*2  l 

input  to  LED/LD  5 

T5 

1  st 

column  of  A  propagates 

VN 

H 

O 

■P 

2nd 

column  of  A  propagates 

to  x4 

3rd 

column  of  A  propagates 

to  T  3 

4th 

column  of  A  propagates 

to  t2 

bi  2 

input  to  LED/LD  1 

b  i  2 

input  to  LED/LD  2 

b+22 

input  to  LED/LD  3 

t>2  2 

input  to  LED/LD  4 

Outputs 


The  arrangement  in  Equation  ( 5-16)  can  also  be  used  to 
perform  a  true  matrix  subtraction  without  the  need  to  handle 
negative  numbers.  For  example,  for 

A  -  B C  =  D  (5-17) 

the  equation  is  rewritten  using  the  earlier  ()+  and  ()“ 
notation  as 

D  =  D+  -  D"  =  (A+  -  A")  -  |(BC)+  -  (BC)"| 

+  +  (5-18) 

=  |  A+  +  (BC)‘j  -  | A'  +  (BC)  | 

Realization  of  Equation  (5-17)  in  the  form  of  Equation  (5-18) 
follows  directly  from  the  matrix  partitioning  used  in  Equation 
(5-16).  The  bipolar  and  matrix  subtraction  algorithm  now  com¬ 
pletes  the  list  of  required  operations  for  Kalman  filtering. 

It  is  left  to  the  final  section  of  this  chapter  to  specify 
the  performance  of  this  processor  completely  as  well  as  to 
compare  its  performance  to  a  software  implementation.  Until 
these  steps  are  accomplished,  the  worth  of  the  processor 
cannot  be  adequately  judged. 

5.3  Design  and  Analysis  of  the  Optical  Kalman  Filter 


It  is  in  this  section  that  the  actual  processor  design 
and  development  subject  to  all  the  previously  outlined  con¬ 
straints  is  reported.  Analysis  of  the  final  design  is  ac¬ 
complished  with  particular  emphasis  on  satisfaction  of  the 
overall  problem  specifications.  The  analysis  also  includes 
a  comparison  of  the  projected  optical  solution  to  one  which 
has  been  formulated  in  software.  It  is  on  the  basis  of 
these  analyses  that  the  optical  architecture  developed  in 
Section  5.2  is  proposed  as  a  viable,  and  perhaps  preferable, 
alternative  to  other  previously  developed  solutions. 

Prior  to  system  design,  an  analysis  of  the  overall 
algorithm  is  required  to  generate  a  minimum  system  processing 
rate  for  the  OKF.  Figure  36  illustrates  the  master  flow 
diagram  of  the  tracker.  For  the  purposes  of  these  analyses, 
each  cycle  of  the  tracker  is  defined  to  start  at  the  time  of 
R(a,0)  availability  from  the  enhanced  optical  correlator  (SOC) 


aster  Tracker  Flow  Diagram 


shown  as  Block  I  in  Figure  36.  Noting  that  the  pointing  sys¬ 
tem  response  must  occur  prior  to  the  next  sampling  time,  some 
assumption  must  be  made  as  to  the  amount  of  time  this  block 
will  require.  For  the  purposes  of  this  research,  the  time 
required  to  generate  and  implement  the  pointing  system  com¬ 
mands  is  assumed  to  be  15  milliseconds.  This  assumption  is 
based  on  two  factors:  (1)  it  is  anticipated  that  implementa¬ 
tion  of  the  tracking  algorithm  developed  in  Chapters  III  and 
IV  will  require  only  small  excursions  (in  terms  of  micro¬ 
radians)  after  each  propagation  cycle  for  correction  of  the 
previously  applied  control  inputs;  and  (2),  since  the  tem¬ 
plate  for  the  EOC  is  generated  in  parallel  with  the  pointing 
commands,  the  template  is  easily  generated  within  this  time 
frame.  With  these  developments,  the  remaining  software  pro¬ 
cessing  time  must  be  specified  to  generate  the  minimum 
throughput  rate  of  the  OKF. 

The  following  computation  times  (11)  were  thus  used  in 


the  analysis. 

Memory  Storage  or  Retrieval  1  usoc 
Addition  or  Subtraction  2.7  usee 
Multiplication  4.1  usee 
Division  6.6  usee 
A/D  or  D/A  Conversion  50  usee 


The  computational  and  memory  manipulation  times  specified 
are  representative  of  single  precision  processing  times  typi¬ 
cal  of  the  IBM  360  series  and  some  smaller  state-of-the-art 
computers. 

The  formal  .  *  _C(a,B)  in  Block  II,  under  these  defini¬ 

tions,  requires  ,  usee.  This  figure  was  obtained  by  64 
iterations  of  th-.  ^  .utine  and  one  iteration  of  the  T2  rou¬ 
tine  defined  in  VII.  For  operations  in  parallel,  the 

determining  time  factor  is  specified  as  the  largest  processing 
time  applicable.  For  example,  since  the  correlation  A/D 
conversion  occurs  in  parallel  with  the  rest  of  the  processing 
steps  of  T1 ,  this  subblock  requires  one  iteration  of  79.7  usee 


and  63  iterations  of  50  usee.  The  other  processing  blocks 
were  analyzed  in  a  similar  manner.  Thus,  the  OKF  measure¬ 
ment  update  cycle  (Block  III  in  Figure  36)  requires  422  ysec 
of  software  processing  time  as  indicated  by  the  routine  out¬ 
lined  in  Table  VIII  (where  N  is  specified  in  Chapter  IV). 

Table  IX  develops  the  processing  time  required  both  on  an 
initial  basis  (formation  of  all  values)  and  thereafter  (forma¬ 
tion  of  only  non-constant  values)  for  the  relinearization 
process  (Block  IV).  These  processing  times  are  613.9  usee 
and  298.9  usee,  respectively.  Finally,  Table  X  lists  the 
task  breakdown  of  the  OKF  propagation  cycle  (Block  V)  which 
specifies  the  486.4  usee  software  processing  time.  Com¬ 
bining  these  requirements  with  the  15  msec  pointing  system/ 
template  formation  (Blocks  VI  and  VII)  time  assumed  earlier, 
the  total  amount  of  non-OKF  processing  time  is  19.7682  msec 
for  the  initial  iteration  and  19.4532  msec  thereafter. 

These  tables  also  specify  the  number  of  optical  pro¬ 
cessing  operations  required  for  the  formation  of  £+ ,  P+, 
and  P~  in  the  OKF.  Adding  these  requirements  from  both 
filter  cycles,  a  total  number  of  194  optical  processing 
operations  will  completely  perform  the  required  linear 
algebra.  By  subtracting  the  non-OKF  processing  time  from 
the  time  per  FLIR  frame  (33.3333  msec),  the  optical  opera¬ 
tions  must  be  accomplished  within  13.5651  msec  on  the  initial 
cycle  and  13.8801  msec  on  subsequent  cycles.  Thus,  the 
minimum  OKF  processing  rate  is  approximately  14.3  kHz. 

With  this  number  in  hand,  the  component  choice  for  the 
optimum  (in  terms  of  both  component  integration  and  uniform 
output  across  the  entire  parallel  signal  representation) 

OKF  can  now  be  specified.  Since  all  of  the  operations  in 
the  OKF  update  cycle  are  bipolar,  the  optical  processor 
must  have  (as  outlined  in  Section  5.2)  23  LEDs/LDs  and  16 
detectors  for  N=M=8.  This  immediately  presents  a  tradeoff 
between  the  LED/LD  power  with  its  resulting  spatial  re¬ 
quirement  and  the  amount  of  acoustic  attenuation  over  the 
resulting  required  (by  the  LED/LD  array)  AOM  aperture.  This 
in  itself  is  the  major  design  constraint. 
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Table  VII.  COM  Determination  Processing  Requirements 


Time 

Operation 

Duration 
( usee ) 

T1 

1-In  parallel,  get  vertical  pixel  number 
get  horizontal  pixel  number,  and  A/D 

9 

convert  one  of  64  pixels 

50 

2-Threshold  pixel  intensity 

2.7 

3-Tempstore  pixel  intensity 

1 

4-Add  value  to  vertical  sum 

2.7 

5-Store  vertical  sum 

1 

6-Retrieve  value  from  tempstore 

1 

7-Multiply  by  vertical  pixel  number 

4.1 

8-Add  to  weighted  vertical  sum 

2.7 

9-Store  in  weighted  vertical  sum 

1 

10-Retrieve  value  from  tempstore 

1 

11- Repeat  4  and  5  for  horizontal  sum 

12- Repeat  6  through  9  for  weighted 

3.7 

horizontal  sum 

8.8 

Total  79.7 

-  =  - 

T2 

1 -Divide  weighted  horizontal  sum  by 

horizontal  sum 

6 . 6 

2-Store  in  memory 

1 

3-Get  weighted  vertical  sum 

1 

4-Repeat  1  and  2  for  vertical  case 

7.6 

Total  16.2 

SIS 


Table  VIII.  OKF  Update  Cycle  Processing  Requirements 


Time 

Operation 

Software 

Duration 

T1 

T 

D/A  convert  H,  P  ,  H  ,  R,  x  , 
and  z  as  per  pipeline 
(one  iteration  of  22  memory 
retrievals  plus  D/A  conversion 
and  seven  iterations  of  D/A 
conversion/memory  retrieval  in 
parallel) 

422  usee 

T2 

Form  H  P"HT  +  R 
(Bipolar  3-matrix  multiply 
with  real-time  resistive 
add,  N=8) 

T3 

Form  HP- 

(Bipolar  2-matrix  multiply, 

N=8) 

T4 

Form  z  -  Hx"  using  Eqn.  (5-18) 
(Bipolar  "2-matrix  multiply 
with  real-time  resistive 
add,  N=8 

T5 

Invert  H P'HT +  R 
(Bipolar  inversion,  N=2) 

T6 

Form  K 

(Bipolar  3-matrix  multiply, 

N  =  8) 

T7 

Form  P+  using  Eqn.  (5-18) 
(Bipolar  3-matrix  multiply 
with  real-time  resistive 
add,  N=8) 

T8 

Form  x+ 

(Bipolar  2-matrix  multiply 
with  real  time  resistive 
add,  N=8) 

Optical 


31  (4N-1) 

23  (3N-1) 

23  (3N-1) 
9  ( 5N-1 ) 

31  (4N-1) 

31  (4N-1) 

23  (3N-1) 


Time 


Operation 


Initial 


Thereafter 


T1 

Zero  $p,  matrices 

(64  elements  each) 

128  ysec 

T2 

Produce  _£p,£p£,  matrices 

(Initial  cycle-60  multiplies, 
11  additions,  77  memory  manip¬ 
ulations,  2  divisions,  and 

2  exponentials  §  50  sec  per 

exponential;  thereafter 

47  multiplies,  10  additions. 

66  memory  manipulations,  and 

2  divisions) 

485.9  qsec 

298.9  ^sec 

Totals 

613.9  u sec 

298.9  usee 

Table  X.  OKF -Propagation  Cycle  Processing  Requirements 


Time 

Operation 

Software 

Duration 

; 

Optical  J 
Operations 

T1 

Form  x"  using  Eqn.  (4-15) 

(8  multiplies  &  8  additions) 

54.4  usee 

! 

T2 

D/A  Convert  $p,  P  ,  and  $p 

as  per  pipeline 
(one  iteration  of  32  memory 
retrievals  plus  D/A  conver¬ 
sion  and  seven  iterations  of 
memory  retrieval/conversion 
in  parallel) 

432  usee 

1 

1 

♦ 

1 

1 

T3 

Form  P" 

(3-maTrix  multiply  with 
real  time  resistive  add, 

N  =8 ) 

2  3  (3N-D 

This  major  design  constraint  defines  the  approach  to 
the  component  selection  stage  of  the  design  process.  Re¬ 
calling  that  the  modulator  is  the  key  component  of  the 
architecture,  it  must  be  specified  first  and  the  other  com 
ponents  (since  there  is  a  greater  range  of  product  avail¬ 
ability  for  each  type)  fitted  to  the  modulator.  However, 
in  order  to  specify  the  modulator,  some  assumptions  must 
be  made  about  the  power  range  (in  terms  of  peak  power) 
required  of  the  source  array.  This  will  allow  an  initial 
design  which  can  be  modified  (if  required)  to  correspond 
to  actual  component  selection. 

From  industrial  data,  the  appropriate  parameters  of 
interest  for  various  AOM  materials  are  listed  in  Table  XI. 

Table  XI.  AOM  Parameters  of  Interest 


Material/  Acoustic  Velocity (Acoustic  Attenuation ' Fig .  of  Merit 
Mode  v  (mm/usec)  .  a0  (dB/us-GH^)  j  M 


PbMO  4/ 
L (001 ) 

TeO  2  / 

L ( 001 ) 

GaP/ 

L  ( 1 1  0) 


0.22 


0.23 


0.8473 


Assuming  for  the  moment  that  each  laser  diode  (due  to  the 
power  versus  size  relationship  and  output  power  required) 
must  be  spaced  1.5mm  center- to-center ,  the  resulting  workin 
aperture  of  34.5mm  when  combined  with  the  AOM  acoustic  velo 
cities  yields  the  following  data  rates  into  the  cell  (acous 
tic  velocity  divided  by  the  laser  diode  spacing)  and  time 
apertures  (working  aperture  divided  by  acoustic  velocity)  . 


Material 


PbMO  1 


Data  Rate  Into  Cell 


2.42  MHz 


2.8  MHz 
4.213  MHz 


’ime  Aperture 


9.5  usee 

8.21  usee 


o.4o  usee 


As  none  of  these  required  data  rates  are  unreasonable  for  the 
associated  electronics,  further  investigation  is  warranted. 
Since  AOMs  are  specified  at  3dB  frequency  response  bandwidths 
(for  the  acoustic  frequencies),  the  requirement  that  only  the 
flat  (i.e.,  uniform)  part  of  the  response  curve  for  a  given 
frequency  range  specifies  that  the  device  response  bandwidth 
be  larger  than  required  just  by  the  necessity  of  producing 
1 6  products  (from  2M  detectors,  M=8)  spaced  at  some  constant 
distance.  As  this  (as  well  as  diffraction  efficiency)  is  a 
function  of  transducer  impedance  matching  (higher  frequencies 
require  smaller  area  transducers) ,  the  additional  desire  is 
that  the  center  frequency  be  kept  as  low  as  possible.  These 

factors  merge  into  specifications  of  10  MHz  frequency  separ¬ 

ation  per  row  product  (150  MHz  total  bandwidth  for  16  pro¬ 
ducts)  ,  a  250  MHz  AOM  acoustic  bandwidth,  and  an  AOM  center 
% 

frequency  of  380  MHz  to  preclude  second  harmonic  effects 
within  the  range  of  cell  frequencies. 

With  these  criteria  defined,  the  total  cell  attenuation 
across  the  AOM  working  aperture  becomes 

Material  a 

PbMO  4  7.55  dB 

Te02  7.47  dB 

GaP  3.00  dB 

where  the  total  cell  attenuation  is  the  acoustic  attenuation 
identified  in  Table  XI  multiplied  by  the  AOM ' s  time  aperture 
and  acoustic  center  frequency.  Combining  this  data  with  the 
superior  thermal  conductivity  of  Gap  (37  times  that  of  Te02). 
the  choice  of  material  is  obvious.  Thus,  a  Gap  AOM  with  a 
center  frequency  of  380  MHz,  a  device  acoustic  bandwidth  of 
250  MHz  (380  MHz  ±  125  MHz),  and  a  transducer  with  appropriat 
shape  to  produce  the  collimated  acoustic  beam  illustrated  in 
Figure  37  is  selected. 


ACOUSTIC  PROPAGATION 


Figure  37  -  Acoustic  Beam  Propagation  in  GaP 

The  entire  OKF  can  now  be  specified.  For  ease  of  imple¬ 
mentation,  the  input  data  rate  is  chosen  to  be  4  MHz  so  that 
custom  clocks  do  not  need  to  be  incorporated  into  the  archi¬ 
tecture.  This  produces  the  specifications  (based  on  A=0.85um 
so  as  to  provide  the  greatest  possible  range  of  laser  diode 
product  selection)  summarized  in  Table  XII,  where  all  values 
are  as  defined  in  this  chapter.  These  specifications  thus 

Table  XII.  AOM  Specifications 


Parameter 

- r 

Value 

Working  Aperture 

LD  Center-to-Center  Spacing 

Attenuation  over  working  aperture 

Bragg  Angle 

Angular  Separation  Between  Products  (A0) 
Undif fracted-Dif fracted  Angular  Separation 
Diffraction  Efficiency  (1  Watt  Drive) 

36 . 34mm 

1  . 58mm 

3.16  dB 

0.442° 

0.023° 

0.712° 

50% 
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require:  (1)  the  detector  selected  must  have  a  minimum  dyna 

mic  range  of  30  dB  (from  the  original  problem  specification) 
at  the  wavelength  of  interest;  (2)  the  laser  diode  chosen 
must  have  a  minimum  dynamic  range  of  33.16  dB ;  and  (3).  the 
lens  system  used  must  have  an  effective  focal  length  of  5m 
[from  Equation  (5-7)3  to  produce  a  center- to-center  spacing 
of  2mm  (due  to  the  minimal  30  dB  dynamic  range  requirement 
plus  an  allowance  for  greater  dynamic  range)  and  a  focussed 
spot  of  145  um  at  the  output  plane. 

Taking  these  criteria  one  at  a  time,  the  detector  cho¬ 
sen  has  the  following  specifications  (1983  Laser  Focus 
Buyer's  Guide,  page  52): 


Minimum  Detectable  Power 
Saturation  Power 
Rise  Time 
Diameter 


500  PW 
1  m  W 
50  nsec 
1.85  mm 


Although  these  detectors  must  be- housed  in  a  custom-built 
array,  their  compact  size  allows  ready  implementation  in 
this  manner. 

Similarly,  the  same  source  (page  184)  defines  the  fol¬ 
lowing  specifications  for  the  laser  diode  selected: 


Peak  Power 

Dynamic  Range 

Repetition  Rate  (Maximum) 

Pulse  Length 

Beam  Divergence 

Wavelength 

Diameter 


20  mW 
54  dB 
5  MHz 
60  nsec 
Collimated 
0.85  um 
1.35  mm 


Again,  the  geometry  is  such  that  a  custom  array  is  easily 
formed. 

This  leaves  only  the  lens  system  to  be  specified. 
Obviously,  a  five-meter  focal  length  lens  resulting  in  a 
10-meter  separation  between  the  A0M  and  the  detection  plane 


would  make  the  system  unwieldy  in  a  physical  sense.  Thus, 
the  sole  method  of  accomplishing  the  required  focal  length 
within  physical  realizable  limits  necessitates  combining  a 
negative  lens  in  tandem  with  a  positive  lens.  This  process 
increases  the  angular  divergence  of  the  products  in  a  linear 
manner,  and  is  most  easily  calculated  in  diopters.  From 
the  1983  Melles  Griot  Optics  Guide,  a  diopter  is  defined 
as  1000  divided  by  the  lens  focal  length  in  mm.  Since  a 
five-meter  focal  length  is  desired,  the  number  of  diopters 
required  is  0.2.  To  keep  the  focal  lengths  of  both  lenses 
at  a  reasonable  limit,  an  arbitrary  boundary  of  T  'Miim  is 
selected.  For  a  negative  lens  with  a  focal  length  of  -400mm, 
the  diopter  conversion  works  out  to  be  -2.5.  This  then 
requires  a  2.7  diopter  positive  lens,  which  works  out  to  a 
370mm  focal  length  lens  by  the  formula  given  above.  Thus, 
all  specifications  are  satisfied  for  the  lens  system,  and 
this  completes  the  specification  of  0KF  componenets. 

To  analyze  the  0KF  completely,  since  it  is  obvious 
from  the  dynamic  range  of  the  laser  diode  that  a  uniform 
response  is  achievable  (across  the  output  parallel  repre¬ 
sentation)  by  introducing  additional  gain  where  required  at 
the  input  (Section  5.2),  it  is  only  necessary  to  calculate 
the  minimum  and  maximum  signals  through  the  system  and  com¬ 
pute  their  resultant  dynamic  range.  Thus,  for  minimum  power, 
the  calculation  is  (using  the  minimum  value  for  each  data 
point) 

(80nW  input)  (.01  diffraction  ef  f  iciency)  (  .  95  lens  trans)=760  p’i 

which  is  above  the  defined  threshold  of  the  detector.  The 
maximum  power  level  through  the  system  is 

(20  mW  input) (.5  diffraction  eff iciency )(. 95  lens  trans)=9.5  mW 

Since  this  value  is  above  the  detector’s  saturation  level,  the 
system  is  specified  to  operate  at  a  maximum  power  level  of 
1  mV/  throughput.  This  yields  an  effective  dynamic  range  of 
61.2  dB  (equivalent  to  20-bit  digital  accuracy)  . 


With  the  entire  OXF  designed  and  analyzed,  the  operation 
of  the  tracking  algorithm  warrants  reinvestigation.  Since  194 
optical  operations  are  required  for  the  filter  linear  algebra, 
the  OKF  accomplishes  all  its  processing  in  48.5  usee  at  the 
4  MHz  data  rate.  This  can  be  contrasted  with  a  software  com¬ 
putation  time  requirement  of  435.2  usee  merely  to  multiply 
a  square  matrix  of  dimension  eight  by  a  vector  of  appropriate 
length .  This  decrease  in  computation  time  results  in  the 
overall  computation  time  of  the  tracker  being  capable  of 
supporting  much  higher  bandwidth  operations. 

For  example,  a  100  Hz  system  bandwidth  could  be  supported 
as  follows.  Given  the  fact  that  the  dynamics  model  used  is 
based  on  a  constant  turn-rate,  it  can  be  postulated  that  for 
a  sample  period  less  than  one-third  the  earlier  definition, 
the  resultant  filter  errors  will  be  reduced  by  about  the 
ratio  since  is  directly  proportional  to  At.  Thus,  by 

merely  maintaining  the  pointing  system  response  rate,  the 
pointing  system  response  time  can  be  cut  to  5  msec  from  the 
earlier  assumption  of  15  msec.  Template  formation  can  also 
be  accomplished  within  this  time  frame.  By  redefining  the 
pointing  system  response  time  in  this  manner,  the  system 
processing  requirements  (optical  and  software)  now  become 
9.8167  msec  for  the  first  iteration  and  9.5017  msec  thereafter. 

This,  however,  would  require  that  the  FLIR  sensor  be 
discarded  in  favor  of  a  custom  array  of  real  time  sensors, 
which  would  be  sampled  every  10  msec.  This  proposal  .nay  pro¬ 
duce  some  very  useful  advantages.  Since  the  sensors  are 
real-time,  smearing  of  the  image  (and  thus  introducing  noise 
into  the  transform)  due  to  relative  motion  between  the  target 
and  tracking  window  trajectories  would  cease  to  exist. 
Additionally,  since  the  noise  mechanisms  of  the  detector 
types  are  different,  the  individual  pixel  noise  term  (here- 
tofor  neglected  due  to  predominance  of  spatially  correlated 
noise)  can  be  substantially  reduced  by  careful  detector  choice. 
As  an  array  of  truly  independent  detectors,  moreover,  the 
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inevitable  crosstalk  between  pixels  is  entirely  negated.  Fi 
nally,  due  to  the  substantially  higher  sampling  rate,  filter 
transient  response  tines  should  decrease  dramatically. 

With  the  entire  tracker  now  completely  specified,  con¬ 
clusions  and  recommendations  can  be  drawn  about  its  worth. 
Thus,  the  next  chapter  relates  the  thesis  research  to  the 
original  goals  expressed  in  Chapter  I  and  projects  these 
relationships  forward  to  their  next  logical  step. 
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VI.  Conclusions  and  Recommendations 
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6.1  Overview 

As  the  major  goal  of  this  research  is  to  demonstrate  that 
the  use  of  optics  and  optical  processing  can  produce  substan¬ 
tial  system  enhancement  at  the  current  state  of  optical 
development,  this  theme  is  central  to  all  conclusions  drawn. 
Coupled  to  this  theme  is  the  necessity  of  identifying  a 
satisfactory  stochastic  algorithm  that  will  be  fully  capable 
of  providing  the  level  of  performance  required  of  a  high 
energy  laser  pointing  system.  These  two  factors  thus  con¬ 
stitute  the  major  bases  for  this  study. 

All  conclusions  and  recommendations  developed  within  this 
chapter  are  therefore  related  to  these  bases.  Although  the 
chapter  is  not  strictly  separated  into  stochastic  and  optical 
sections,  care  is  taken  to  identify  each  conclusion  and  recom¬ 
mendation  made  with  respect  to  its  relationship  to  both  the 
stochastic  and  optical  bases.  By  examining  all  factors  in 
this  light,  all  research  reported  in  this  thesis  is  tied 
together  to  produce  an  end  product  which  is  usable  in  either 
its  current  form  or  easily  modified  to  reflect  alternate 
performance  criteria. 

Section  6.2  thus  identifies  all  appropriate  conclusions 
which  can  be  drawn  from  the  research  performed.  Recommenda¬ 
tions  developed  from  these  conclusions  are  then  listed  in 
the  final  section  of  this  chapter.  In  this  manner,  the  thesis 
is  concluded. 

6.2  Conclusions 

Since  two  trackers  using  the  same  stochastic  algorithm 
have  been  proposed  in  this  thesis,  both  the  FLIR-constrained 
and  unconstrained  trackers  must  be  examined  in  turn.  In  this 
manner,  any  differences  in  conclusions  can  be  brought  to 
light,  and  will  allow  recommendations  to  be'  developed  based 
on  these  differences. 
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Initially  examining  the  FLIR-constrained  tracker,  it  is 
important  to  note  that  the  use  of  optical  processing  in  this 
tracker  produces  no  substantial  performance  enhancement. 

This  is  chiefly  due  to  the  fact  that  the  optical  Kalman  fil¬ 
ter  and  enhanced  optical  correlator  have  system  processing 
bandwidths  far  in  excess  of  that  provided  by  the  use  of  the 
FLIR  array  sensor.  This  limitation  thus  results  in  the 
tracking  system  being  required  to  idle  a  minimum  of  40* 
between  measurements  and  represents  a  substantial  waste  of 
the  capability  of  the  tracker's  individual  components.  On 
a  stochastic  basis,  given  the  previous  research  (8,21,23), 
the  developed  algorithm  can  be  reasonable  anticipated  to 
yield  acceptable  performance  against  the  postulated  targets 
of  interest.  Therefore,  should  the  FLIR-constrained  tracker 
be  the  solution  of  choice,  it  must  be  concluded  that  a  soft¬ 
ware  implementation  of  the  developed  stochastic  algorithm  is 
more  than  adequate  and  indeed  preferable  to  process  the  data. 

However,  if  the  unconstrained  tracker  is  used  to  develop 
an  alternative  solution,  the  use  of  the  optical  processing 
techniques  described  within  this  thesis  become  mandatory  to 
provide  the  necessary  system  bandwidths  to  support  the  over¬ 
all  tracker  throughput  rate.  Additionally,  use  of  the  same 
stochastic  algorithm  at  the  higher  measurement  rate  can 
reasonably  be  expected  to  produce  enhanced  performance  in  the 
critical  areas  of  transient  response  times  and  error  statis¬ 
tics,  as  outlined  in  Chapter  V.  Given  the  statistics  developed 
in  the  earlier  research  (8,21,23),  it  is  indeed  possible  that 
the  unconstrained  tracker  can  produce  a  pointer-constrained 
system  due  to  its  higher  processing  bandwidth  against  most 
targets  of  interest.  Thus,  it  must  be  concluded  that  the 
unconstrained  tracker  developed  within  this  thesis  is  indeed 
the  preferable  solution. 

Finally,  some  conclusions  which  are  common  to  both  of  the 
trackers  proposed  in  this  thesis  must  be  identified.  Given 
the  state-of-the-art  of  integrated  optics  (27),  either  tracker's 
optical  Kalman  filter  could  be  implemented  in  this  form.  This 

1 1 1 


could  produce  major  savings  in  size,  weight,  and  power  consume 
tion  requirements  for  the  system.  Additionally ,  the  use  of 
optics  in  either  case  does  permit  the  use  of  more  sophisticate 
stochastic  algorithms,  particularly  if  a  multiple-model  adap¬ 
tive  filter  is  to  be  used.  While  these  conclusions  do  not 
seem  to  be  critical  for  the  case  of  a  ground-based  laser  wea¬ 
pon,  translation  of  the  laser  system  to  either  an  airborne  or 
space  environment  may  very  well  require  a  combination  of 
these  benefits,  chiefly  provided  by  the  use  of  optical  pro¬ 
cessing  techniques. 

In  sum,  it  is  the  conclusion  of  this  thesis  that  the 
unconstrained  tracker  developed  within  this  thesis  is  the 
solution  of  choice  for  a  ground-based  laser  weapons  system. 
Since  any  successful  implementation  of  this  tracker  requires 
the  use  of  the  optical  processing  techniques  developed  and 
analyzed  within  this  thesis  within  the  constraint  of  using 
currently  available  optical  components,  the  major  basis  of 
the  thesis  has  been  satisfied.  Lastly,  given  the  other  con¬ 
clusions  developed  within  this  section,  the  use  of  the  iden¬ 
tified  optical  processing  techniques  will  provide  substan¬ 
tially  more  flexibility  in  the  choice  of  stochastic  algorithms 
to  match  the  scenario  of  application.  As  the  use  of  sto¬ 
chastic  algorithms  is  indicated  due  to  their  enhanced  per¬ 
formance  over  correlation  algorithms,  the  processing  bed 
developed  is  applicable  to  the  widest  range  of  problem 
scenarios . 


6.3  Recommendations 

The  conclusions  developed  in  Section  6.2  yield  three 
recommendations  for  future  research  efforts.  These  recommen¬ 
dations  are  based  not  only  on  the  conclusions  but  also  on 
the  premise  that  the  maximum  benefit  be  derived  from  the 
work  performed  within  this  thesis. 

The  first  recommendation  is  that  a  complete  Monte  Carlo 
performance  analysis  be  performed  for  both  trackers  proposed. 
This  performance  analysis  should  be  conducted  in  two  chases. 


First,  both  trackers  should  be  evaluated  against  the  trajec¬ 
tories  outlined  in  the  truth  model  presented  in  Chapter  II. 

After  careful  evaluation  of  the  tuning  parameters  involved, 
both  trackers  should  then  be  analyzed  using  the  same  techniques 
against  trajectory  tapes  of  actual  aircraft  and  missile  per¬ 
formance.  In  this  manner,  realistic  statistics  can  be  genera¬ 
ted  while  gaining  tuning  experience  with  the  filter. 

Secondly,  a  test  bed  optical  Kalman  filter  should  be 
developed  and  used  to  generate  a  true  modulation  transfer 
function  after  all  optical  components  have  been  optimized  as 
identified  in  Chapter  V.  This  test  bed  can'  then  be  used 
to  evaluate  the  system  throughput  rate  of  any  stochastic 
algorithm.  Upon  combining  the  results  of  this  recommendation 
with  those  of  recommendation  one,  a  reasonable  projection  of 
filter  performance  versus  sampling  rate  variations  can  be 
developed  for  a  variety  of  different  algorithms. 

Lastly,  it  is  strongly  recommended  that  a  serious  inves¬ 
tigation  into  alternatives  to  the  FLIR  array  sensor  be  initiated 
As  the  major  system  limitation,’  the  FLIR  array  sensor  is  the 
major  stumbling  block  to  enhanced  performance.  Development 
of  a  reasonable  alternative  to  the  FLIR  sensor,  particularly 
in  time  response  characteristics,  will  then  enable  the  develop¬ 
ment  of  a  final  optimal  package. 

It  is  believed  that  the  results  of  these  recommendations 
will  enable  the  final  specification  of  a  pointing  and  tracking 
system  which  will  yield  enough  performance  benefit  to  warrant 
implementation.  It  is  also  postulated  that  the  experience 
gained  from  both  this  thesis  and  implementation  of  its  recom¬ 
mendations  will  allow  a  much  easier  translation  of  the  final 
result  from  a  ground-based  version  to  other  scenarios.  In 
sum,  the  optimization  of  the  final  product  is  the  end  goal  of 
this  thesis,  its  conclusions,  and  its  recommendations. 
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A ocendix  A 


Two-Dimensional  Finite  Discrete  Fourier  Transform 
Interpretation  and  Test 


Similar  to  the  presentation  of  J.  W.  Goodman,  Reference 
6,  the  Rect  function  is  defined  to  extract  one  period  of  the 
spatial  intensity  function  as 


RectN(x,y) 


1  0<  x  <  N-1  ,  0<  y  <  N-1 

( A  —  1  ) 

0  otherwise 


Reproducing  Equations  (2-3)  and  (2-4)  of  Reference  23  and 
using  the  Rect  function  to  define  the  area  sequence  g(x,y) 
as  being  zero  outside  the  interval  0  £  x  £  N-1 

g ( x ,  y )  =  g  '  ( x ,  y )  Rectjy  ( x ,  y )  (A-2) 


G(f  ,f J 


x  y  Lx=0  y=0 

g(x,y) 


N-1  N-1  :4rIL(f  x  +  f  y)‘ 

I  Z  g  ( x ,  y )  e  N  x  y 


RectN(fx,f  )  ( A—  3 ) 


N-!  N-1  ^(V  +  V*’ 

Z  £  G  ( f  ,  f  )  e  w  x  y 

f  =0  f  =0  x  y 
w  x  y 


Rect,T(f  ,f  ) 
N  x  y 


( A  —  4) 


G(f  ,f  )  is  the  Finite  Discrete  Fourier  Transform  of  g(x,y). 

*’'•  y 

The  Rect  function  of  Equation  (a-1)  is  separable  in  the  two 


independent  variables: 

Rectjj  ( f  x  >  f  y )  -  RectN(fx)  Rect^(fy) 
Equation  (A-3)  can  now  be  written  as 


( A-  5 ) 


0(W 


n-i 


l  X(f  ,y)  exp 

Ly=o  x 


Rect,,  ( f  ) 

«  y 


( A  -  6 ) 


where 


x(fx»y)  = 


N-1  ^(O)' 

Z  g(x,y)  exp 
x  =  0 


ReotN(fx) 


( A-  7) 


Equation  (A-7) ,  X(f  ,y)  corresponds  to  an  N-point  one  dimen¬ 
sional  Discrete  Fourier  Transform  for  each  value  of  the  row 
index  y.  X(f^,y)  is  the  result  of  N  one  dimensional 


transforms,  one  for  each  row  of  g(x,y).  In  figure  A-la,  if 
y  is  held  constant,  say  y=2,  and  a  one-dimensional  Fourier 
Transform  is  accomplished  across  that  row,  the  variation  in 
the  image  intensities  from  column  to  column  would  result  in 
the  zero  frequency  of  that  variation  being  located  at  x'=1 
and  y'=2  of  Figure  A-lb,  while  the  coefficient  associated 
with  the  fundamental,  the  first  harmonic,  in  both  directions 
would  be  found  at  x’=2  and  y,=2  with  its  conjugate  located 
at  x ' =N  and  y'=2. 


(a)  Original  Data  Array 


DC  T, 


h - 

2)  T  (2 

•  N) 

A 

✓ 

AV 


(b)  Transformed  Data 

Figure  A-1  -  Row  Transform  Result  of  2-D  DFT 


118 


To  complete  the  two-dimensional  Discrete  Fourier  Trans¬ 
form,  Equation  (A-6)  expresses  how  to  implement  the  iJ  one¬ 
dimensional  transforms  along  each  column.  To  summarize  this 
general  interpretation  of  the  two-dimensional  Discrete  Fourier 
Transform,  Equations  (A-6)  and  ( A - 7 )  show  how  the  two-dimen¬ 
sional  transform  is  achieved  by  using  a  one-dimensicnal  trans¬ 
form  on  the  rows  first  and  then  on  the  columns,  or  vice  versa. 

For  the  case  where  the  image  intensity  function  is  sep¬ 
arable  in  the  two  independent  variables  x  and  y,  for  example 
g(x,y)  having  the  property  that 

g(x,y)  =  gi(x)  g2(y)  ( A- 8) 


the  two-dimensional  Discrete  Fourier  Transform,  G(f  ,f  ), 

X  y 

becomes  the  product  of  the  one-dimensional  independent  trans¬ 


forms  Gi(f  )  and  G2(f  ). 


G(f  ,f  )  =  Gi(f%)  G2 ( f  ) 
x  y  x  y 


( A  -  9 ) 


A  function  of  two  independent  variables  is  separable  within  a 
specific  coordinate  system  if  it  can  be  written  as  a  product 
of  two  functions  of  one  independent  variable  each.  The  two- 
dimensional  transform  degenerates  to  holding  the  row  index 
constant,  for  example,  and  running  the  one-dimensional  trans¬ 
form  across  the  columns,  for  any  one  of  the  N  rows.  Similarly 
the  column  index  is  held  while  the  one-dimensional  transform 
is  accomplished  across  the  rows,  and  this  is  done  for  any  of 
the  columns.  The  results  of  these  two  independent  transforms 
are  then  multiplied  together. 

This  leads  to  some  interesting  algebra  which  can  be  ex¬ 
ploited  to  test  the  implementation  of  the  two-dimensional 
Discrete  Fourier  Transform.  In  Figure  A-1,  T^  2)  con¬ 

sist  of  the  product  of  the  two  one-dimensional  transform 
fundamental  coefficients  for  that  row  or  column.  If  x  is 
the  column  coordinate  and  y  is  the  row  coordinate,  then 
T,  \  is  defined  as 


(2 


fist  harmonic  in  y 

1  =v 

1 st  harmonic  in  x 

!  for  column  2 

-yO,2) 

• 

for  row  2 

r  X 


(1,2 


( 2  ,N ) 


(N-1, 2) 


(N-1 ,N) 


1st  harmonic  i 
for  column 


in  yl  - 

mn  N  j=y(1  ,N 


conjugate  oi  1st 
harmonic  in  x 
for  row  2 


'conjugate  of  1st' 
harmonic  in  y 
for  column  2 

'conjugate  of  1st' 
harmonic  in  y 
for  column  N 


iy(1 ,2) 


:y(1,N) 


1st  harmonic  in  x  _ 
for  row  N-1 


(A- 11  ) 


conjugate  of  1st' 
harmonic  in  x 
for  row  N-1 


J'A(1 ,N-1 

( A  -  1  2 ) 

EX(1 ,N-1) 
(A-13) 


if  yi  <  N/2,  xi  <  N/2 


A ( y i »xi) 


(yi-1 )  harmonic  in  y 

for  column  x _y  (  y  i  -  1  ,  xi ) 


( xi-1 )  harmonic 
for  row 


lie  in  x 
3w  y  i  J  ' 


( x  i  - 1  ,  y  i ) 


If  the  variation  in  intensity  across  every  row  is  equal  to 
the  variation  in  intensity  for  every  other  row  and  the  vari¬ 
ation  across  all  of  the  columns  is  similarly  set  equal,  then 
the  componenets  of  the  transform  can  be  defined  as 


y(1 ,2)  =  a  + 

y(1,M)  =  a  +  jb 

y(1 ,2)  =  a  -  J'b 

y(1,N)  =  a  - 


(1,2) 


=  c  +  jd 


c  +  jd 


X(1 , N - 1 )  =  c 

x(1,2)  =  jd 


( A- 1 4) 


y(l  ,N)  a  "  J'b  X(  1  ,N- 1 )  =  c"  J’d 
Using  Equation  C A - 1 4 )  to  solve  (A-10)  through  (A-13) 

T(2,2)=y(1,2)x(1,2)  =  ^a+J'b^c  +  ^d^=^ac'bd^+^<'bc  +  ad^  C  A  —  1  5 ) 

T(2,N)=y(1,N)x(1,2)=(a+jb)(c-jd)=(ac+bd)-j(  ad~bc)  ( A  - 1 6 ) 

T(N-1 ,2)  =y(  1  ,2)x(1,N-1)=^a_J'b^c  +  Jd^  =(ac  +  bd)  +J  (ad- be)  (A-1  7) 

T(N-1  ,N)  =y(1  ,N)X(1  , N - 1 )  =  ^a“J'b^  (c- jd)s(ac-bd)  -j (ad+bc)  (A-1  8) 

Equation  (A-1 5 )  for  is  the  conjugate  of  Equation  (A-1 8 ) 

for  T ^ j  under  the  conditions  of  uniform  variations  on 
any  row  and  similarly  for  the  columns,  which  just  implies 
separability . 


■  i."  j  j  ■ 


To  implement  the  two-dimensional  Fourier  Transform,  sub¬ 
routine  Fourt,  a  common  Fortran  Subroutine,  was  used.  To 
test  subroutine  Fourt  to  see  where  it  places  harmonics  with  a 
two-dimensional  array,  the  results  of  Equations  ( A - 1 5 )  through 
(A-18)  were  used.  Figures  A-2  through  A-7  show  the  results 
of  these  tests. 


£F**! 

^  -  ■' 


Figure  A-2  -  Cos  2ttx/6 


1  21 


Figure  A-1  had  as  its  input  a  24x24  array  with  magnitudes 
of  each  element  varying  only  across  the  rows  via  cos(2ttx/6). 
Using  the  above  interpretation  of  a  two-dimensional  Fourier 
Transform,  the  fundamental  frequency  assumed  would  be  one 
which  corresponded  to  one  period  across  the  24-element  array. 
Since  the  input  obviously  fits  four  cycles  in  that  space,  the 
only  nonzero  component  of  the  Fourier  Transform  expected  would 
be  at  the  fourth  harmonic.  The  transform  would  also  be  expec¬ 
ted  to  be  real  only  since  the  input  is  a  pure  cosine.  Figure 
A-2  is  then  encouraging  in  that  the  only  nonzero  components 
are  where  expected  and  the  imaginary  portions  of  that  answer 
are  extremely  small.  The  reason  that  there  are  nonzero  ele¬ 
ments  only  in  the  first  row  is  that  when  the  transform  was 
accomplished  across  the  columns,  there  was  no  variation  in  the 
y-coordinate ,  which  is  equivalent  to  taking  a  transform  of  a 
constant  one.  That  transform  results  in  a  one  in  the  zero 
frequency  components  of  each  column  which  is  the  first  row  and 
zero  everyplace  else.  When  the  multiplication  of  Equation 
(A-9)  is  accomplished.  Figure  A-2  results. 

Figure  A-3  had  an  input  of  (cos  2irx/6)(cos  2my/24)  which 
is  the  same  variation  in  the  x-direction  and  a  variation  in 
the  y-direction  which  corresponds  to  exactly  one  period.  The 
only  nonzero  result  from  the  transform  along  the  columns,  vari¬ 
ation  in  y,  should  be  along  the  fundamental,  which  it  is. 

Figures  A-4  through  A-7  are  just  further  examples  to 
demonstrate  Fourt.  Figure  A-4  had  an  input  of  the  fourth  har¬ 
monic  in  the  x-direction  and  the  fundamental  in  the  y-direc¬ 
tion.  Figure  A-5  contains  the  sixth  harmonic  in  the  x-direc¬ 
tion.  Figure  A-6  contains  the  twelfth  harmonic  in  the  x-di¬ 
rection  which  only  appears  as  a  conjugate.  Figure  A-7  con¬ 
tains  only  the  twelfth  harmonic  in  both  directions  and  en¬ 
closes  spatial  frequencies  within  boxes.  A  thorough  under¬ 
standing  of  these  results  is  needed  to  understand  the  pro¬ 
cesses  of  taking  derivatives  and  shifting  in  the  transform 
domain . 
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In  summary,  Figure  A-8  is  provided  to  show  the  output  of 
subroutine  Fourt.  The  DC  component  is  the  product  of  the  zero 
frequency  components  of  the  one-dimensional  transforms  in  both 
directions.  Elements  ^  through  24)  can  ke  ^-nter~ 

preted  similarly  to  Equations  (A-10)  through  ( A— 1 3) • 
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Figure  A-8  -  Output  of  Fourt 


jendix  B 


Derivation  of  the  Matrix 


The  derivation  of  the  linearized  F  and  matrices  for 
the  nonlinear  constant  turn-rate  target  acceleration  dynamics 
model  is  not  as  simple  as  it  is  for  a  linear  filter.  This  is 
because  the  point  of  linearization  is  dependent  on  the  current 
value  of  the  state  vector  and  cannot  be  precomputed.  Repeat¬ 
ing  Equation  (4-16) 


F(t^)  =  3_f  ["x,  tj  /3x 


x=x(ti  ) 


(B-1) 


and  recalling  the  definition  of  the  nonlinear  function 
f|x(t),t|  from  Equation  (4-12),  the  F  matrix  has  the  form 


F(t.)  = 


( B-2) 


where 


FI  =  -(oj/A1)(A2  -  4wx^  +  2x^x^,) 

F2  =  2x3w(x5A1  +  2x4A2)/A1z 
F  3  =  2x3x4w/A1 
F4  =  -2x3u)/A1 

F5  =  -2x4oj(x6A1  +  2x4A2)/A12 
F6  =  -u>j>  -  2x4(x^A1  +  2x4A2)/A12J 
F7  =  2x2uj/A1 
F8  =  -2x-,x,w/A1 


in  which 


A1  =  x,  +  xz, 


A2  =  -  x,x 


u)  =  A2/A1 


=  correlation  time  of  the  atmospheric  jitter 

and  with  the  components  of  the  state  vector  specified  in 
Chapter  Iv  written  as 

T 


X  s  [x1  X 


2  x3  X4  X5  x* 


x7  x8^ 


Once  the  F(t^)  matrix  is  determined,  the  state  transition 
matrix  can  be  evaluated  using  the  quasi-static  approximation 
of  Equation  (4-19)  except  for  the  two  atmospheric  terms 
which  are  obtained  in  exact  closed  form.  Specifically,  this 
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Kalman  filter.  The  dynamics  model,  based  on  a  constant 
turn-rate  target  acceleration  model,  is  deemed  to  be  a  better 
representation  of  the  true  target  dynamics  than  a  linear 
first-order  Gauss-Markov  target  acceleration  model  for  the 
cases  of  interest.  Optical  processing  techniques  are  com¬ 
pletely  specified  for  the  correlation  stage  to  perform  the 
required  correlations  in  real  time,  and  the  filter  stage  to 
perform  the  linear  algebra  required  for  the  Kalman  filter. 
This  extensive  use  of  optics  allow  the  development  of  two 
tracking  algorithms  based  on  the  same  models:  a  FLIR-con- 
strained  tracker  with  a  30  Hz  frame  rate  and  an  unconstrained 
tracker  with  a  100  Hz  frame  rate  using  real-time  sensors  in 
place  of  the  FLIR. 
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